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ABSTRACT 


The  inviscid  barotropic  vorticity  equation  is  integrated  under  a 
variety  of  assumed  initial  and  boundary  conditions  corresponding  to  linear 
and  nonlinear  box  modes,  forced  nonlinear  box  modes,  north  wall  forced 
modes,  and  linear  and  nonlinear  Rossby  waves  (with  and  without  mean 
advection).  The  former  two  classes  of  problems  are  defined  within  a 
closed  domain.  The  latter  two  are  partially  or  totally  open  to  a specified 
external  environment  and  therefore  represent  prototype  limited-area 
calculations  for  the  ocean. 

To  determine  the  extent  to  which  the  accuracy  and  efficiency  of 
limited-area  calculations  depend  on  the  numerical  integration  scheme, 
each  test  problem  is  solved  independently  using  the  finite -difference  (FD), 
finite -element  (FE)  and  pseudospectral  (PS)  techniques.  The  three 
numerical  models  differ  primarily  in  the  formal  accuracy  of  their  spatial 
approximations  and  their  treatment  of  vorticity  at  outflow  points  along  the 
boundary.  The  FD  model  employs  a centered  second- order  differencing 
scheme  and  requires  an  extrapolatory  (computational)  boundary  condition 
to  fix  the  values  of  vorticity  at  outflow  boundary  points.  The  FE  model, 
which  represents  i|r  and  ^ as  a summation  of  piecewise  linear  elements, 
is  of  fourth  order  for  the  linearized  one- dimensional  advective  equation. 
Further,  a technique  is  developed  by  which  the  determination  of  the 
interior  values  of  Q is  decoupled  from  that  of  the  boundary  values;  hence, 
the  vorticity  boundary  conditions  can  be  implemented  without  iterative 
techniques.  Lastly,  the  "infinite  - order  " PS  model  avoids  the  assumption 
of  lateral  periodicity  by  expanding  i|r  and  Q in  a double  series  of 
Chebyshev  polynomials.  The  resulting  vorticity  equation  is  solved  in 
the  spectral  domain  using  a modified  alternating  direction  implicit  method. 

All  three  models  are  of  second  order  in  time  and  have  conservative 
formulations  of  the  nonlinear  terms. 

Integrations  of  moderate  length  (5-10  periods  of  the  known  analytic 
solution)  are  performed  to  determine  the  accuracy,  stability  and  efficiency 
of  each  model  as  a function  of  problem  class  and  the  associated  physical 
and  computational  nondimensional  parameters.  The  most  important  of 
these  parameters  are  e,  the  Rossby  number,  v,  the  number  of  spatial 
degrees  of  freedom  (grid  points,  expansion  functions,  etc.  ) per  half  wave- 
length of  the  reference  solution,  and  q,  the  number  of  time  steps  per 
period  of  the  reference  solution.  The  latter  two  parameters  are  non- 
dimensional  measures  of  the  spatial  and  temporal  resolution  of  the 
numerical  approximation. 

These  tests  show  that  all  three  models  are,  in  general,  capable  of 
delivering  stable  and  efficient  solutions  to  linear  and  weakly  nonlinear 
problems  in  open  domains  (Os  es  0.4,  4s  vS  10,  64s  qs  128).  Despite 
their  added  complexity,  however,  the  FE  and  PS  models  are  on  the  average, 

4 and  15  times  more  accurate  respectively  than  the  FD  model  even  taking 
into  account  its  increased  efficiency.  The  results  also  suggest  that  given 
a judicious  selection  of  a frictional  (filtering)  mechanism  and/or  computational 
boundary  condition  (to  suppress  the  accumulation  of  gridscale  features), 
each  of  the  models  can  be  made  similarly  accurate  for  highly  nonlinear 
calculations  (e  » 0.4). 
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1.  INTRODUCTION 

Limited-area  or  regional  modelling  of  the  oceans  is  important  for  a 
number  of  scientific  and  potentially  practical  problems.  These  include 
real-time  forecasting/hindcasting  studies  of  the  oceanic  interior  and 
regions  of  intense  near-surface  currents  (such  as  the  Gulf  Stream),  and 
idealized  local  dynamic  studies  of  those  subregions  of  the  oceanic  gyre 
thought  to  contribute  to  the  generation  and  maintenance  of  the  mean  and 
transient  fields  of  motion.  The  former  studies  relate  to  ongoing  and 
anticipated  measurements  from  a variety  of  modern  techniques  including 
satellite  surveys  of  the  ocean  [1,2].  The  latter  are  prompted  by  recent 
advances  in  our  modelling  and  understanding  of  the  mesoscale  features  of 
the  ocean  circulation  [3]. 

Numerical  eddy- resolving  general  circulation  models  in  gyre-scale 
systems  (EGCM's)  have  shown  that  transient  eddies  similar  to  those 
observed  play  an  important  role  in  determining  the  mean  flow  [4-6].  The 
very  high  horizontal  resolution  needed  to  resolve  adequately  the  dynamics 
on  the  scale  of  the  eddies,  however,  make  EGCM's  expensive  to  run  even 
when,  in  practice,  basins  of  only  a few  eddy  wavelengths  in  horizontal 
extent  are  used.  Furthermore,  gyre -scale  models  are  complicated  by  the 
fact  that  they  contain  many  subregions  of  distinct  physics  [7,  8], 

For  statistically  homogeneous  subregions  of  the  gyre,  local  dynamical 
studies  can  be  made  with  periodic  boundary  condition  models.  Such  models 
assume  that  the  physics  is  locally  determined  and  essentially  independent 
of  information  such  as  scales  and  amplitudes  which  could  be  generated 
elsewhere  and  continuously  transported  across  the  boundaries.  Such 
"process"  models  have  been  used  to  investigate  the  dynamical  properties 
of  the  mesoscale  eddy  field  under  simulated  mid-ocean  conditions  in 
regions  well  removed  from  boundary  layer  effects  [9,  1 0 ].  Many  subregions 
of  the  gyre  obviously  violate  the  above  assumptions.  For  such  regions,  as 
well  as  for  many  other  regional  hydrodynamic  problems  other  more 
complicated  boundary  conditions  are  necessary  [H].  Such  problems 
include  limited-area  oceanic  forecasting  [12],  coastal  modelling  [13],  and 
the  study  of  intense  current  systems  and  their  associated  instabilities. 
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The  determination  of  valid  and  convenient  forms  of  boundary 
conditions,  particularly  at  points  of  outflow,  constitutes  a major, 
essentially  unresolved,  problem  in  the  modelling  of  many  hydrodynamic 
systems  over  regional  domains.  The  choice  of  boundary  conditions 
involves  a number  of  physical,  mathematical,  and  numerical  (or  computa- 
tional) considerations.  On  large  and  regional  scales,  the  conditions  should 
correctly  represent  or  parameterize  the  interaction  of  the  (arbitrary) 
volume  of  fluid  with  its  surroundings.  Smaller  scale  physical  phenomenon 
generated  internally  within  the  region  should  not  be  trapped  but  allowed  to 
pass  out  of  the  domain;  i.  e.  , the  model  boundary  should  be  transparent 
for  such  small-scale  processes.  The  mathematical  problem  consisting  of 
model  equations  and  boundary  conditions  should  be  well  posed.  The 
numerical  scheme  chosen  for  computational  purposes  should  be  of  desired 
accuracy  and  acceptable  efficiency. 

In  its  computational  algorithms,  a given  numerical  scheme  may 
involve  the  use  of  boundary  information  which  is  not  in  principle  required 
in  the  well-posed  analytical  problem.  Such  auxiliary  information  is  known 
as  a computational  boundary  condition  and  should  be  chosen  for  convenience 
and  efficiency,  but  should  not  in  principle  affect  physical  results.  In 
practice,  specific  problems  appear  to  require  special  conditions  and  a 
trial  and  error  approach  is  usually  indicated.  Alternate  choices  of 
boundary  conditions  do  effect  the  accuracy  and  stability  of  regional 
computations.  For  instance,  in  viscous  limited-area  calculations  it  is 
known  that  incorrect  specification  of  boundary  conditions  or  data  generally 
leads  to  the  evolution  of  boundary  layers  adjacent  to  outflow  [14,  15], 

Several  types  of  boundary  conditions  have  been  used  to  avoid  the  problems 
associated  with  outflow,  including  extrapolatory  formulas  [16,17],  and 
radiation  conditions  [11],  In  inviscid  systems,  Charney,  Fjortoft  and  von 
Neumann  [18]  originally  argued  heuristically  that,  given  the  stream- 
function  at  all  boundary  points,  only  the  values  of  vorticity  at  inflow  points 
were  needed  to  determine  consistently  solutions  to  the  barotropic  vorticity 
equation.  This  conjecture  has  since  been  proven  more  rigorously  [19]. 

In  practice,  however,  unless  iterative  or  implicit  numerical  techniques  or 
one-sided  differencing  schemes  are  used,  inviscid  calculations  require 
some  auxiliary  relationship  by  which  to  prescribe  vorticity  at  outflow 
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boundary  points.  Such  a numerical  scheme  is  often  sensitive  to  the 
specific  choice  of  this  computational  boundary  condition.  An  example  of 
this  has  been  given  by  Shapiro  and  O'Brien  [20]  in  the  context  of  a limited- 
area  meteorological  forecast.  They  showed  that  the  method  of  character- 
istics worked  well  as  a computational  boundary  condition  but  that  the 
specification  of  vorticity  at  outflow  boundary  points  on  the  basis  of  known 
boundary  data  led  to  numerical  instability. 

In  this  report,  we  evaluate  the  feasibility  of  performing  barotropic 
open  ocean  calculations  by  investigating  and  comparing  the  accuracy, 
efficiency  and  stability  of  three  limited-area  numerical  models  based 
respectively  on  the  finite -difference,  finite -element  and  spectral 
approximation  methods.  The  physical  boundary  conditions  used  are  the 
Charney-Fjortoft-von  Neumann  conditions.  The  calculations  are  mostly 
inviscid  but  in  some  cases  a dissipative  filter  is  included.  The  codes 
differ  primarily  in  the  details  and  accuracy  of  their  spatial  discretization 
schemes  and  in  their  treatments  of  the  vorticity  at  outflow.  The  statement 
of  the  nondimensional  vorticity  equation  and  a description  of  the  numerical 
techniques  are  given  in  Sections  2 and  3. 

The  three  models  have  been  tested  and  intercompared  for  a variety 
of  prototype  physical  problems  in  closed  and  open  basins  and  over  a range 
of  the  nondimensional  physical  and  computational  parameters  corresponding 
to  each  problem  class.  First,  the  unforced  (homogeneous)  solutions  to  the 
linear  and  nonlinear  vorticity  equation  in  a closed- basin  are  found  and 
compared  to  the  known  exact  and  perturbation  solutions  for  linear  and 
nonlinear  box  modes  respectively  (Section  4).  With  the  addition  of  a body 
force,  various  exact  nonlinear  closed-basin  solutions  are  constructed  and 
tested  (Section  5).  Next,  oscillations  driven  by  time -dependent  distribution 
of  boundary  values  along  the  northern  (open)  boundary  are  examined 
(North  wall  forced  modes  - Section  6).  Lastly,  in  a fully  open  domain, 
linear  and  nonlinear  Rossby  wave  solutions  with  and  without  mean  ad- 
vection  are  obtained  and  intercompared  (Section  7).  Model-model  inter- 
comparisons of  this  type  have  been  carried  out  for  simple  advective 
problems  in  closed  basins  (e.  g.  , Orszag  and  Israeli  [21]).  To  our 
knowledge,  however,  this  is  the  first  such  study  that  encompasses 
limited-area  hydrodynamic  modelling  problems  as  well. 
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The  development  of  a regional  qoasigeostrophic  modelling  capability 
has  two  immediate  intended  applications.  Firstly,  forecasting/hind- 
casting studies  of  the  mid-ocean  mesoscale  eddy  field  will  be  undertaken. 
For  this  work  the  MODE  and  POLYMODE  synoptic  data  sets,  can  be 
utilized  in  a variety  of  ways  both  as  boundary  information  and  forecast 
verification  data.  Secondly,  motivated  by  recent  EGCM  simulations  of 
eddy  generation  and  eddy  mean  field  interaction  in  regions  of  intense 
mean  flow,  the  finite -amplitude  instabilities  of  idealized  jets  will  be 
investigated  as  a function  of  environmental  parameters  and  inlet  conditions. 
Ultimately,  if  possible,  real  Gulf  Stream  meander  region  observational 
data  will  be  used.  For  these  purposes,  baroclinic  models  will  be  needed. 
The  extension  of  the  following  results  to  a multi-level  configuration  is, 
however,  straightforward.  A two-level  version  of  the  code  described  in 
Section  3.  3 is  already  operational. 
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2.  MODELLING  EQUATIONS,  METHODOLOGY  AND 
FORMAT  OF  RESULTS 

We  consider  the  barotropic  vorticity  equation  on  a (3-plane  which 
can  be  written  in  dimensional  form  as 

a ) 0 £ x £ L 

+ J(*,  ) (C  + f)  = F(x,y)  x (1) 

> OSysL 

y 


where 


and 


f = fn  + Py 


j(*  ) = J>i  ±. 

lT»  ’ 9x  9y 


In  addition, 


C = v* 


it  JL 

9y  9x 


(2) 


is  the  relationship  between  streamfunction  (f)  and  vorticity  (£),  and 
F(x,  y)  represents  the  effects  of  a body  force,  if  any.  Dissipation  has  been 
neglected  for  three  reasons.  First,  the  inviscid  system  poses  a simpler 
physical  and  numerical  problem,  one  for  which  many  analytic  and  per- 
turbation solutions  are  available.  This  is  the  basis  of  our  testing  of  the 
limited-area  models  described  in  Section  3.  Second,  quadratic  conservation 
laws  are  available  for  non- dissipative  physical  and  numerical  systems. 

This  property  also  contributes  to  the  evaluation  of  model  performance. 
Lastly,  by  ignoring  explicit  higher -order  friction,  we  sidestep  for  the 
moment  the  question  of  the  correct  specification  of  vorticity  boundary 
conditions  on  outflow,  which  are  not  formally  needed  for  integration  of  the 
inviscid  system.  The  assumption  of  inviscid  dynamics  does,  however, 
require  that  greater  care  be  taken  to  construct  a numerical  scheme  which 
is  stable  in  the  absence  of  explicit  dissipative  (that  is,  smoothing) 
mechanisms.  As  we  shall  see,  such  filtering  is  in  fact  necessary  to  main- 
tain stability  in  some  cases. 
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If  we  now  nondimensionalize  x,  y,  t,  and  i|i  by  d,  d,  (pd)  and 
(Vgd)  respectively,  then  (1)  becomes,  in  nondimensional  form, 


+ eJ(+,  ) V2i|i  + tx  = F(x,y) 


0 s x s x. 


OS  ys  y. 


where  the  Rossby  number 
e = VQ/pd2 

and 

XB  = Vd  yB  = Ly/d  • 

The  parameters  d and  Vq  are  taken  to  be  the  characteristic  length  and 

velocity  scales  of  the  anticipated  field  of  motion.  Note  that  the  length 

scale  d does  not  correspond  to  the  basin  dimensions  L or  L ; hence, 

x y 

Xg  and  yB  are,  in  general,  greater  than  one. 

The  modelling  strategy  developed  herein  involves  the  integration  of 
Eq.  (3)  for  several  sets  of  initial  and  boundary  conditions  corresponding 
to  succeedingly  more  complex  physical  phenomena  in  closed  and  open 
domains.  The  problems  we  will  consider  include  linear  and  nonlinear  box 
modes,  forced  nonlinear  box  modes,  linear  and  nonlinear  north  wall  forced 
modes  (meander-induced  forcing),  and  linear  and  nonlinear  Rossby  waves. 

The  sequence  of  linear  problems  (box  modes,  north  wall  forced  modes, 
and  Rossby  waves)  serve  as  pivotal  calculations  for  which  no  boundary 
values  of  vorticity  are  formally  required.  With  the  addition  of  nonlinearity, 
both  accuracy  and  stability  of  model  calculations  can  be  assessed  as  a 
function  of  € for  closed-domain  problems  in  which  strict  conservation 
laws  apply  (nonlinear  and  forced  nonlinear  box  modes),  and  partially  open 
and  totally  open  domain  problems  in  which  interaction  with  the  surrounding 
environment  is  possible  and  the  question  of  computational  boundary 
conditions  arises  (nonlinear  north  wall  forced  modes  and  nonlinear  Rossby 
waves).  The  former  experiments  are  the  most  easily  understood.  The 
latter  series  of  tests  - particularly  the  nonlinear  Rossby  waves  with  mean 
advection  - are  those  most  relevant  to  future  open  ocean  modelling 
applications. 

1 1 
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Each  problem  thus  defined  has  been  solved  using  three  independent 
(and  quite  different)  numerical  techniques.  The  finite -difference  (FD) 
scheme  is  second-order  accurate  in  space  and  time  and  has  the  advantage 
of  being  easily  coded.  It  should,  in  addition,  be  the  most  efficient  of  the 
three  models  for  a fixed  number  of  spatial  degrees  of  freedom.  The  finite- 
element  (FE)  method,  though  somewhat  more  complicated  than  the  FD 
scheme,  is  known  to  be  of  fourth  order  for  certain  advective  problems  [22]. 
In  general,  we  expect  smaller  errors  with  the  FE  model,  but  at  a slightly 
higher  computational  cost.  Lastly,  pseudospectral  (PS)  approximation 
techniques  [21]  offer  greatly  reduced  spatial  truncation  errors  in  comparison 
to  both  of  the  other  two  methods.  The  PS  model  is,  therefore,  formally  the 
most  accurate  but  thereby  it  may  be  subject  to  instabilities  not  seen  in  the 
FD  formulation  - see,  for  instance,  Section  7.  Of  the  three  models,  it  is 
also  the  most  difficult  to  code  (although  it  can  be  made  comparably  efficient 
if  care  is  taken  to  optimize  the  spectral  transforms). 

Since  analytic  or  perturbation  solutions  are  available  for  many  of  the 
prototype  physical  problems  examined  herein,  direct  measures  of  numerical 
error  are  available  for  each  model.  Of  particular  interest  are  the  RMS 
errors  in  streamfunction  and  vorticity,  and  the  normalized  difference  in 
integrated  kinetic  energy  as  a function  of  time;  these  area  integrated  error 
measures  are  defined  as 


RMS(( ')  =j 

ff»c  - MJ^a'2  -^j172 

(4a) 

RMS(C')  = | 

JK  - «a’2  <MK>2  -I'72 

(4b) 

and 

NDIF(NRG)  = 

where  subscripts  c and  a refer  to  the  computed  and  analytic  solutions 
respectively  and  a primed  quantity  represents  a difference  from  the 
reference  solution.  Using  error  measures  of  this  sort,  it  is  possible  to 
ascertain  the  accuracy  of  each  model. 


{}}|v*c|2dA-  JJl?tJ2dA}/jjK 


dA 


(4c) 


Table  I summarizes  the  results  of  all  the  experiments  as  a function 
of  problem  type  and  the  associated  nondimensional  parameters.  The  first 
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seven  columns  of  Table  I refer  to  the  experiment  number  and  the  (not 
necessarily  independent)  quantities 


(i) 

e 

Rossby  number 

(ii) 

XB  = yB 

nondimen sional  basin  size 

(iii) 

1 

N 

number  of  spatial  degrees  of  freedom  in  each 
direction 

(iv) 

h = xb/(N-1) 

nondimensional  mesh  interval 

(v) 

A 

number  of  half  wavelengths  or  turning  points  of 
the  reference  solution  within  the  domain  (non- 
dimensional  measure  of  the  structure  of  the 
solution) 

(vi) 

v = (N-l)/A 

nondimensional  spatial  resolution  (number  of 
degrees  of  freedom  per  turning  point  of  the 
reference  solution) 

(vii) 

At 

nondimensional  time  increment 

(viii) 

T) 

nondimensional  temporal  resolution  (number  of 
time  steps  per  period  of  the  reference  solution). 

The  last  four  columns  tabulate  the  duration  of  the  experiment  (in  periods) 
and  the  final  values  of  the  three  error  measures  defined  above.  The 
duration  of  simulations  which  suffered  numerical  instability  are  denoted  by 
brackets.  No  RMS  error  values  are  listed  for  these  experiments.  Inter- 
mediate columns  of  Table  I are  reserved  for  special  parameters  representative 
of  each  problem  class.  These  will  be  introduced  in  Sections  4 through  7. 

Accompanying  Table  I are  a series  of  figures  which  show  in  more  detail 
the  results  of  one  experiment  for  each  problem  category.  Figures  (1-3)  are 
typical.  The  first  two  figures  give  contour  plots  of  Q and  and  Q'  and 
\|r 1 respectively  at  the  end  of  the  simulations  for  each  model.  Figure  3 shows 
the  corresponding  variation  with  time  of  RMS(C')  and  RMS(f')  for  each 
model.  Contouring  intervals  and  scaling  information  are  given  in  the  figure 
captions.  Note  in  particular  that  the  RMS  error  curves  are  not  necessarily 
scaled  similarly,  even  within  a given  figure. 

The  former  (contour ) plots  give  a visual  estimate  of  the  wavenumber  con- 
tent of  the  computed  solution  (andhence  an  estimate  of  its  formal  spatial  accuracy) 
and  reveal  any  localized  features  (physical  or  numerical)  that  might  occur. 

This  is  particularly  important  in  those  nonlinear  experiments  for  which  only 
approximate  or  linearized  reference  solutions  are  available  and  for  which  the 
RMS  errors  are  therefore  difficult  to  interpret. 
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3.  MODEL  FORMULATION 
3.  1 Finite-difference  model 


A traditional  discrete  formulation  of  Eq.  (3)  results  by  approximating 
all  derivatives  by  second-order  centered  finite  differences.  The  advantages 
of  finite  differences  are  the  simplicity  of  the  required  coding  and  our 
accumulated  experience  of  using  this  classical  technique.  With  centered 
(leapfrog)  time  differencing,  the  vorticity  and  Poisson  equations  become: 


Q 


k+l 


£ 


k-l 


J*(*k,Ck) 

h 


At 

h 


(5a) 


and 


6 f 

xx 


k+l 


6 A 

yy 


k+l 


h2Ck+1 


(5b) 


. k k k k 

where  finite -difference  operators  6 (i|f  ) = i|r.  , . -f.  , . and  6 (^  ) = 

k k k ^ ^ j J ^ j J xx 

i|r. , i . - 2ij(.  . + f . . . The  Arakawa  [23]  Jacobian  J*  is  given  by  the 

finite- difference  approximation  to  the  equivalent  form  1/3 f [VXC  - t £x]  + 

+ t(Ky)x  - (Kx)yJ  + [-(tyS)x  + (♦  xC)y]}  • This  expression  conserves 

vorticity,  energy  and  enstrophy  when  integrated  over  a closed  domain. 

The  Poisson  equation  (5b)  with  Dirichlet  boundary  conditions  is  approximated 

by  the  standard  five  point  Laplacian  operator  and  solved  with  the  NCAR 

cyclic  reduction  subroutine  package.  Finite  differences  are  of  second-order 

2 2 

accuracy,  so  the  total  discretization  error  is  0(h  + At  ). 

For  linear  problems  (e  = 0),  vorticity  on  the  boundary  E does  not 
enter  the  problem  in  either  the  vorticity  equation  or  the  Poisson  equation. 

For  nonlinear  problems  (e  / 0),  vorticity  is  specified  at  inflow  points 
according  to  the  Charney-Fjortoft-von  Neumann  boundary  condition  [18]. 
Centered  finite  differences  require  boundary  data  everywhere,  and  in 
contrast  to  the  analytic  problem,  some  auxiliary  relationship  (that  is,  a 
computational  boundary  condition)  must  be  assumed  at  outflow  points  in 
order  to  determine  the  vorticity  there  (unless  an  iterative  technique  is  used 
to  fix  £ on  £). 


Optimally,  the  computational  boundaries  of  an  open  ocean  model 
should  be  transparent  to  signals  impinging  on  them.  Thus,  the  formation 
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of  boundary  layers  on  or  wave  reflections  from  the  boundary  are 
undesirable.  The  most  successful  computational  boundary  condition 
considered  in  this  study  is 


c 


k 

B 


k+1 

B- 


+ c 


k-1 

B-l 


(6) 


which  was  introduced  by  Sundstrom  [16].  Here  B,  B-l,  B-2  represent 
a boundary  point  and  its  1st  and  2nd  normal  interior  neighbors.  Davies 
[17]  demonstrated  the  stability  properties  of  this  closure  for  a variety  of 
nonlinear  problems. 

There  are  several  possible  physical  interpretations  of  statement 
(6):  either  (a)  it  is  equivalent  to  equating  the  time  and  spatial  averages 
of  Q at  point  B-l  (a  kind  of  smoothness  condition  at  outflow); 

(b)  Eq.  (6)  is  equivalent  to  = c £ where  c = Ax/ At  (a  "local" 
wave  equation);  or  (c)  it  is  a low- order  spatial  extrapolation  scheme. 

In  order  to  implement  the  Sundstrom/ Da  vies  formula,  the  quantity 

is  eliminated  by  application  of  the  vorticity  equation  evaluated  at 
the  point  B-l,  t=  kAt.  This  yields  an  implicit  set  of  equations  for  the 
boundary  vorticity.  Formally,  this  requires  the  inversion  of  a hepta- 
diagonal  matrix  with  cyclic  ordering  of  the  points.  However,  it  can  be 
shown  that  by  elimination  of  those  boundary  points  which  are  corner 
neighbors,  a simple  tridiagonal  system  results,  provided  there  is  at 
least  one  inflow  point  (Appendix  I).  Note  that  the  boundary  vorticity  is 
calculated  after  the  interior  vorticity  and  the  streamfunction. 

Other  computational  boundary  conditions  investigated  in  this  study 

k k k 

are  the  Kreiss  extrapolation  £g  = 2£^  ^ - £^  2 [1 1 ] and  the  condition 


3.  2 Finite-element  model 

In  the  finite -element  formulation,  we  assume  a set  of  basic  functions 
consisting  of  piecewise  polynomials,  the  simplest  being  piecewise  linear 
elements  arranged  in  a rectangular  lattice.  In  one  dimension,  each 


A 


r 
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element  is  a chapeau  or  hill  function,  and  in  two  dimensions,  a pyramidal 
function,  each  centered  on  the  lattice  point  with  base  width  2h  [22], 

Then  the  basis  functions  have  the  property 


0jzj 

1 if  p = q 

= ' 

(7) 

q -p 

0 if  p / q 

where  z is  the  p 
-P 

-th  lattice  point.  All  fields  can 

be  expressed  in  terms 

of  a summation  of  the  basis  functions;  for  instance 

} 

♦ (x,  t) 

= E +q(t)  V?) 

(8) 

C(x,t) 

= E 

(9) 

where  x is  a general  point  within  the  domain. 

The  weights  ijr ^ and  are  obtained  by  the  Galerkin  procedure. 
Substitute  (8)  and  (9)  into  the  vorticity  equation,  multiply  by  a basis 
function  and  integrate  over  the  entire  region.  Since  each  element  over- 
laps those  adjacent  to  it,  in  general  each  equation  will  contain  terms 
from  its  eight  neighbors. 

For  closed-domain  problems,  we  use  centered  time  differences. 
The  resulting  finite -element  form  of  the  vorticity  equation  is 

M(Ck+1)  = M(Ck-1)  - 2 AtQk  = Pk+1  (10a) 

where 

Qk  = \ J*(*k’£k)  - Ih  w(y>6xHk)  • 

h 


In  partially  and  completely  open  domains,  however,  a second-order 
Adams -Bashforth  time -differencing  scheme  is  used  to  avoid  computational 
instability.  In  such  cases, 
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The  mass  matrix  M has  the  entries  M - ff 0 0 dA.  M is  factorable 

(x)  (v)  * ^ P ^ 

into  two  parts,  M = W ® W 7 , where  ® denotes  tensor  multiplication 

(x)  ( V ) 

and  the  matrices  W and  W'7  are  tridiagonal  with  the  local  form 
1/6(1  4 1).  The  operator  ® can  be  interpreted  to  mean 

successive  multiplications- -first  rowwise,  then  columnwise --by  the 
matrix  W = = W^.  Note  that  the  superscripts  (x)  and  (y)  refer 

to  the  order  in  which  the  tridiagonal  multiplications  are  done.  At  each 
time  step,  therefore,  Eq.  (10)  can  be  written 

M(ck+1)  = w(x)  ® w(y)(Ck+1)  = Pk+1 

k+1 

where  Q is  the  N x N matrix  of  values  of  Q at  time  step  (k+1).  In 
this  form,  it  is  clear  that  Eq.  (10a)  is  equivalent  to  two  tridiagonal  matrix 
systems  each  of  size  N x N.  The  mass  matrix  can  therefore  be  readily 
inverted.  Note  that  if  W is  set  equal  to  the  identity  matrix,  Eq.  (10a) 
reduces  to  the  finite-difference  form  (5a).  The  Jacobian  term  is  precisely 
the  Arakawa  Jacobian  employed  with  finite  differences.  In  fact,  the 
Arakawa  form  is  derivable  from  the  finite -element  formulation  [24]. 

Fix  [22]  has  shown  that  linear  elements  for  the  linearized  advective 
equation  + U£x  = 0 produce  fourth-order  accurate  phase  errors.  To 
maintain  this  accuracy  for  the  vorticity  equation,  the  solution  of  the 
Poisson  equation  for  the  streamfunction  must  also  be  of  fourth  order. 

This  is  accomplished  by  the  method  of  deferred  corrections  [25].  Note  that 

K(i|r ) = hVfr  + J2  h4|?4  - 2 2 ) » + 0(h6) 


where  K is  the  usual  five- point  Laplacian 


K 


Therefore  two  successive  Poisson  solutions  yield  f to  fourth  order  inthe 

following  manner.  First,  obtain  a second- order  solution,  ijr.,  from 

2 1 
K(|j)  = h £ . Then,  a fourth-order  estimate  of  | is  the  solution  to 
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K(t)  = 


h2r  4.  h_ 

h C + i2 


V2C  - 2 


2 2 
3x  dy 


(11) 


Finite  elements  require  more  computational  work  per  time  step 
than  do  finite  differences.  First,  two  tridiagonal  inversions  must  be 
performed  to  determine  the  vorticity  field,  contrasted  with  a simple 
direct  substitution  in  finite  differences.  Second,  two  calls  must  be  made 
to  the  Poisson  solver  instead  of  one.  The  significant  increase  in  accuracy 
plus  the  virtue  of  using  a technique  based  on  a variational  principle  justify 
this  increased  computational  effort  in  many  applications,  as  indeed  they 
will  here  (see,  for  instance,  Section  7). 

In  the  finite -element  model,  vorticity  boundary  conditions  are 
implemented  in  the  following  manner.  For  ease  of  presentation,  we 
introduce  three  types  of  points  and  their  respective  computational  mole- 
cules m,  that  is,  their  local  contribution  to  mass  matrix  M: 


(a)  interior  points 


mj  = 1/36 


(12a) 


(b)  regular  boundary  points  m^,  = 1/36 


and 


(c)  corner  points 


mSE  = 1/36 


/ 0 0 0 \ 

\ 0 2 4 /. 


(Eastern  wall) 
(12b) 


(Southeast 

corner) 


(12c) 


The  lattice  point  associated  with  the  given  element  is  denoted  by  the  under- 
line. Analogous  computational  molecules  exist,  for  regular  boundary 
points  on  the  northern,  southern  and  western  walls  and  for  the  southwest, 
northwest  and  northeast  corner  points. 

Assume  first  that  vorticity  is  specified  everywhere  on  the  boundary 
(corresponding  to  inflow  everywhere)  and  solutions  are  needed  only  for 
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the  interior  points.  It  is  then  easy  to  show  that  this  is  equivalent  to  the 
system  of  equations: 


M4(Ck+1)  = P4 


where 


for  interior  points 


£b/36  at  (eastern)  boundary 


/ 0 0 1\ 

I 0 0 4 ) Cn/36  at  (southeast)  corner 

Vl  4 1/  B 


(13) 


(14) 


is  the  specified  boundary  vorticity,  and  M4  = ® W4y\  where 

W „ (=  = W(y))  is  the  (N-2)  X (N-2)  tridiagonal  matrix 

4 4 4 


1 4 


(15) 


The  subscript  (4)  refers  to  the  corner  terms  in  (15). 

Next  consider  the  case  where  vorticity  is  not  specified  anywhere  on 
the  boundary,  as  in  a basin  totally  enclosed  by  solid  and/or  outflow  sides. 
Here  solutions  are  sought  for  the  entire  field  (interior,  boundary  and 
corner  points).  It  is  then  easy  to  show  that  the  combined  system  including 
contributions  of  the  form  (12a-c)  is  equivalent  to 

M2(Ck+1)  = P2  <16> 


(17) 


is  now  N x N. 
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These  two  cases  represent  extreme  situations  - -either  all  inflow, 
or  all  outflow  (and/or  solid  boundaries).  We  seek  a method  which  will 
allow  a general  mix  of  inflow  and  outflow  sections  around  the  boundary. 
We  can  do  this  by  defining  a new  matrix,  which  has  the  following 

computational  structure 

/ m_  for  points  not  at  or  adjacent  to  the 

j boundary  (i,  j/  1,2,  N-1,N) 


mI  - 2 mE 


1 1 .1 
ml"  2mE  " 2mS+  4mSE 


for  interior  points  adjacent  to  the 
(eastern)  boundary  but  not  near  a 
corner  (i  = N-l.  j / 2 or  N-l) 

for  interior  corner  (southeast)  points 
(i  = N-l,  j = 2) 

(18) 


and  so  on  for  points  adjacent  to  other  boundaries  and  corners. 

It  can  be  shown  that  this  new  formulation  decouples  the  determination  of 
the  interior  vorticity  from  that  of  the  boundary  vorticity  and  is  equivalent 
to 

k+t 

M7/,(C  ) = P-7/o  (19) 


where  P7/7  has  the  same  relationship  to  P as 


to  M given  in 


(18).  In  addition,  M7 /7  = W 


7/2  - w 7/2  ® " 7/2 
7/2  1 


where 


1 7/2 


and  all  of  the  unknowns  are  interior  points.  In  short,  we  use  the  known 
dynamic  relations  between  boundary  and  interior  points  to  disconnect  the 
solution  of  the  one  from  the  other. 

Furthermore,  given  the  interior  values  from  the  inversion  of 
each  of  the  four  boundaries  decouples  from  the  other  boundaries 
and  all  of  the  corners- -for  instance  by  defining  a new  matrix 

ME  = ME  " 2 MSE  “ 2 MNE  1 ^ 
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at  the  eastern  boundary.  The  solution  for  the  boundary  values  requires 
4 inversions  (one  for  each  wall)  of  a tridiagonal  system.  The  corner 
point  values  then  follow  algebraically.  In  the  finite- element  approxi- 
mation, the  order  of  calculation  is  therefore  the  following:  interior 
vorticity,  vorticity  at  regular  boundary  points,  vorticity  at  corner  points, 
and  lastly  the  streamfunction.  The  reader  should  note,  however,  that  this 
solution  procedure  does  not  strictly  guarantee  that  values  of  vorticity  on, 
or  adjacent  to,  points  of  inflow  be  in  exact  dynamic  balance. 

3.  3 Pseudospectral  Model 

We  seek  a (discrete)  spectral  solution  to  Eq.  (3)  subject  to  some 
appropriate  set  of  boundary  and  initial  conditions.  For  definiteness, 
consider  specifying  boundary  values  of  streamfunction  and  vorticity  in 
the  manner  first  suggested  by  Charney,  Fjortoft  and  von  Neumann  [18]. 

That  is,  we  take  as  given  quantities  the  values  of  j everywhere  along 
£,  and  Q at  those  points  along  E characterized  by  mass  influx. 

Boundary  values  of  vorticity  at  outflow  points  are  therefore  unconstrained; 
they  are  computed  as  part  of  the  calculation.  Under  these  boundary 
conditions,  a completely  enclosed  domain  is  a special  case,  one  for 
which--owing  to  the  absence  of  any  inflow  at  all- -vorticity  need  never 
be  specified  at  any  time  along  E.  In  analogy  to  the  analytic  problem,  we 
make  a computational  distinction  between  problems  contained  within 
bounded  regions  and  those  characterized  by  partially  or  fully  open  domains. 

3.  3.  1 Closed  domain 

In  a closed  system,  the  advective  terms  in  Eq.  (3)  are  treated 
explicitly  by  a leapfrog  time -differencing  scheme.  Under  this  second- 
order  approximation,  the  vorticity  and  Poisson  equations  become,  in  the 
usual  notation, 


ck+1  = f-1  . 2ntjcJ(tk,£k)  + ^j 

= Rk+1(x,y) 

(21) 

vV+l  = Ck+1  . 

(22) 

and 
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In  space,  we  adopt  a pseudospectral  approximation  technique  for  which 
the  dependent  variables  are  expanded  in  a series  of  Chebyshev  polynomials} 
that  is,  let: 


and 


N 

N 

tk(x,y)  = 

E 

E 

n=0 

m=0 

N 

N 

Ck(x,y)  = 

E 

E 

n=0 

m=0 

1 

N 

N 

R (X,y)  = 

E 

E 

n=0 

m=0 

re 

A 

X = 

(2x  - 

xfi)/x 

A 

y = 

(2y  - 

yB)/y 

(23a) 


(23b) 


(23c) 


-1 


and  Tn(x)  = cos(n  cos  x)  is  the  Chebyshev  polynomial  (of  the  first 


kind)  of  degree  n,  a function  of  the  linearly  stretched  coordinates  x 


and  y.  On  the  collocation  grid  (x^,  y^)  = (cos(7Tp/N),  cos(7rq/N)), 
Os  p,  qsN,  this  implies  for  instance  that 


*k(W  * J Jo  Vncos  jtr  jc°5  jT* 


) Z2E 


J ffqm 


demonstrating  the  important  fact  that  a Chebyshev  transform  is  equi- 
valent to  a cosine  transform  on  the  (non-uniform)  collocation  grid 


(Xp,  y ) and  as  such  can  be  implemented  very  efficiently  using  special 
forms  of  the  fast  Fourier  transform  algorithm  [26], 


Under  these  definitions,  Eqs.  (21)  and  (22)  can  be  rewritten  in 


bk+1  = Rk+1 
nm  nm 


and 


[a3™  + a^]^1  = bk+1 
1 1 nm 


nm 


j , ai 

nm’ 

1U  XV 

nm 

0 s n, 

m £ N 

(24) 

0 s n, 

m s N 

(25) 

99ft* 
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yy  y v 

where  a and  a'*  satisfy 
nm  nm  7 


N N 


8x2 

la 

n=0 

la 

m-0 

N 

N 

XX 

a 

nm 

E 

E 

n=0 

m-0 

a2 

N 

N 

3 

2 

E 

E 

3y 

o 

a 

m-0 

N 

N 

ayy 

nm 

E 

E 

n=0 

m=0 

T (x)  T (y) 
n'  m 77 


Unlike  periodic  models,  the  resulting  spectral  scheme  can  accommodate 
quite  arbitrary  boundary  conditions  on  E . In  the  present  bounded 
geometry,  i (r(x,  y)  is  given  on  the  boundary  by  some  set  of  values,  let  us 
say  In  terms  of  the  spectral  coefficients  this  is  equivalent  to 

requiring  that 


E E - *E(+1’>rq' 

n=  0 m=  0 n M 


0 s q £ N 


N N 

E (-1)“  E a T (y  ) = 

n=  0 m=  0 M 


♦ 2(- y ) OsqsN 


N N 

E E a t (x  ) 

m=0  n =0  nm  n p J 


■ *z<v+1) 


Os  ps  N 


(26a) 


(26b) 


(26c) 


N 

E <-*> 

m-0 


E a m Tn(x_)  = OspsN 

^0  nm  n p I rE  p7 


(26d) 


These  conditions  are  imposed  on  the  Poisson  equation  (22)  by  using  the 
spectral  analogue  of  the  tau  method  [27],  that  is  by  neglecting  the  highest 
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order  dynamic  equations- -those  for  n = N-l,  N and  m - N-l,  N- -in 
(25).  The  remaining  equations  are  then  supplemented  by  boundary 
conditions  (26a-d),  written  in  their  equivalent  Chebyshev  series  form, 
to  close  the  problem.  The  resulting  matrix  equations  are  not  sparse; 
however,  they  are  quite  easily  diagonalized.  The  details  of  the  solution 
have  been  given  by  Haidvogel  and  Zang  [28]  who  show  that,  for  sufficiently 
smooth  fields  and  given  accuracy,  solutions  to  Poisson's  equation  can  be 
computed  at  least  as  efficiently  by  these  spectral  techniques  as  by  certain 
second  and  fourth-order  finite-difference  methods. 

k+1 

Once  the  spectral  coefficients  a have  been  determined,  thus 
yielding  i|r  (x,  y)  at  the  next  time  level,  the  velocity  components  u =-+  y 
and  v = +^i  can  be  computed  from  well  known  Chebyshev  derivative 
relations.  These  in  turn  are  combined  to  give  the  nonlinear  term 
J(t,  C)  = V • (v£)  by  the  simple  pseudospectral  procedure 


<Yi>lp,q  = W-p,yq)£(ip,y<])]+^M*p,yq)C(ip,yq)) 


where  the  product  vQ  is  determined  locally  by  physical  space  multipli- 
cations on  the  collocation  grid  (x  , y ) but  the  derivatives  d/dx  and 
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8/9y  are  performed  spectrally.  The  resulting  scheme  is  of  infinite- 
order  accuracy  and  can  be  constructed  so  as  to  conserve  any  of  the 
higher  order  invariants  such  as  energy  and  enstrophy;  however,  it 
retains  the  effects  of  high  wavenumber  aliasing  at  full  strength  [29]. 

The  effects  of  aliasing  can  be  identically  removed,  but  at  a large  cost 
in  computational  efficiency  (approximately  a factor  of  2). 


3.3.2  Open  domain 

In  the  spectral  approximation,  boundary  conditions  can  only  be 
correctly  applied  if  the  highest-order  terms  are  treated  in  some  way 
implicitly.  For  those  problems  in  which  the  domain  of  integration  is  not 
bounded  by  an  impermeable  surface  and  for  which,  therefore,  values  of 
the  vorticity  are  to  be  specified  at  points  of  inflow,  this  implies  that  the 
advective  term  must  be  treated  somewhat  differently  than  in  a closed 
region.  To  do  so,  we  adopt  the  following  alternating  direction  implicit 
(ADI)  time -differencing  scheme. 


A 
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^k+1/2  + CM  j__  (UoC)k+l/2  = Li(ck+l/2}  = Ck.i|i^(voal 


{ev  [(v-vQ)C]  + ix}k  + ^ fev-  [(v-v0)C]  + ijix)k  1//2 


(27a) 


k+1  eAt  9 , , vk+1  _ L ,rk+l.  _ -k-1/2  evt  _8_(  c)k+l/2 

£ +—  - L2(C  ' ~ Q 2 8x  lu0U 

- ^ {ev-  [(v- vQ)C ] + +xJk+1//Z  + {ev  [(v-v0)£]  + i|rx3k 


whe 


uk+1(x,y)  = (^){uk+1(l,y)  - uk+1(-l,y)}  + uk+1(-l,y) 
vk+1(x,y)  = (^-)(vk+1(x,  1)  - vk+1(x,  -1)]  + vk+1(x,  -1) 


(27b) 


uk+1(-l,y),  uk+1(  + l,y),  vk+1(x,+l)  and  vk+1(x, -1) 


are  the  known  distributions  of  normal  velocity  at  time  step  (k+1)  along 
the  western,  eastern,  northern,  and  southern  boundaries,  respectively. 
In  effect,  this  semi-implicit  procedure  removes  and  treats  implicitly 
that  portion  of  the  advective  term  which  arises  from  contributions  due  to 
non-zero  normal  velocities  at  the  domain  edges.  Besides  assuring 
stability  of  the  computational  scheme,  the  splitting  of  the  advective  term 
relaxes  the  restrictive  Courant  condition  which  arises  for  explicitly 
differenced  inflow/outflow  problems  due  to  the  crowding  of  Chebyshev 
collocation  points  near  the  domain  boundaries. 

The  solution  of  each  half  step--Eq.  (27a)  or  (27b)- -proceeds 
similarly.  Consider  (27a).  The  implicit  advective  effects  introduce  a 
coupling  only  the  x direction.  In  fact,  along  any  line  y = y = constant, 


u*+1(x,  y ) is  at  most  linear  in  x.  Under  these  circumstances,  operator 

0 ’ q 
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.k+1/2 


) can  be  expressed  spectrally  as  a sequence  of  tridiagonal 


matrix  equations,  one  for  each  y (Os  qs  N).  Similar  remarks  hold  for 
k+1  1 

operator  (C  ) which  is  repeatedly  solved  along  lines  of  constant 
x = Xp  (Os  pS  N).  In  either  case,  vorticity  boundary  conditions  are 
selectively  introduced  in  place  of  higher  order  (dynamical)  equations  in 
each  tridiagonal  system.  Zero,  one  or  two  vorticity  conditions  are 
imposed  depending  on  the  corresponding  number  of  end  points  of  the  line 
x = Xp  or  y = y^  which  are  inflow  points.  The  resulting  numerical 
scheme  fixes  the  value  of  Q only  at  inflow  points;  however,  all  values 
of  boundary  vorticity  are  in  exact  dynamic  balance  with  the  interior. 

(This  is  not  true  of  either  the  FD  or  FE  models--see  the  preceding 
sections.  ) Note  that  problems  involving  inflow  only  in  one  coordinate 
direction- -such  as  the  north  wall  forcing  problem  to  be  discussed 
shortly- -can  be  handled  in  a single  step. 

Once  Eqs.  (27a)  or  (27b)  yield  £^+l/2  or  £^+1^  the  associated 

k+l/2  k+1 

Poisson  problem  for  the  streamfunction  fields  t)  ' or  i|i  are 
solved  as  outlined  above  for  a closed  basin. 
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4.  LINEAR  AND  NONLINEAR  BOX  MODE  TESTS 


In  this  and  the  following  sections,  we  briefly  describe  the 
formulation  and  selected  results  of  the  prototype  numerical  tests 
mentioned  in  Section  1.  For  a more  complete  summary  of  the  results, 
the  reader  is  referred  to  Table  I and  Figures  1-27.  (See  also  Section  2, 
pages  5-8.  ) 

4.  1 Formulation 

The  class  of  exact  solutions  to  the  linear  vorticity  equation  (e  = 0) 
satisfying  homogeneous  streamfunction  boundary  conditions  on  E are  the 
box  modes  or  normal  modes  of  the  basin.  These  can  be  written 


| (x,  y,  t)  = sin(Xx)  sin(py)  cos(x+t/2) 


Os  xS  Xg 

0sysyB 


(28) 


where  x and  t have  been  scaled  with  respect  to  d and  (pd)  , 
respectively  and  d is  taken  to  be  the  scale  length  of  the  travelling  wave 
(wavelength/27T).  The  parameters  X and  p,  and  the  domain  size  Xg  = y^ 
are  related  to  the  integer  mode  numbers,  m and  n,  by  the  relations 

. /#  2 . 2.1/2 

X = m/(m  + n ) 

/,  2 . 2,1/2 
p = n/(m  + n ) 
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and 

T|  = 4ir/&t  = temporal  resolution  parameter. 

As  mentioned  previously,  we  use  N to  refer  to  the  number  of  spatial 
degrees  of  freedom  (grid  points,  spectral  functions,  etc.)  which 
characterizes  the  x and  y discretization  of  each  model. 

For  e « 1,  we  can  obtain  an  approximate  solution  to  the  non- 
linear box  mode  problem  using  a small  amplitude  expansion  in  e.  Let 

V(x,  y,  t)  = t0  + Gilfj  + e2\|f  2 + . . • 

where  (|q  is  the  linear  box  mode  solution  (28).  Then,  to  first  order 

*2*»  + *i*  = -J(Vv2*o)  <29> 


which  has  the  solution 

tj(x,  y,t)  = sin(2py)  j - ^ sin2  \x  + 


| cos(t  + 2x) 

4(l+2|x^)  I 


sinh  R(xg-x) 
sinh  Hxg 


cos(t  + x/2) 


(30) 


with 

R = (4 p2  - l/4)^2,  p > 1/4  . 

Taking  the  expansion  to  second  order,  the  right  hand  side  of  the 

equation  for  if  ^ has  a component  proportional  to  This  secularity 

destroys  the  uniform  convergence  of  the  approximation  for  large  t. 

Following  Pedlosky  [30],  it  can  be  shown  that  by  introducing  the  new  time 
2 

scale  t = t(l  + e 6),  all  forcing  terms  proportional  to  can  be 

suppressed  for  a suitable  choice  of  6.  The  perturbation  solution 

tfr  = ^o  + €^l  can  ^lerefc>re  be  corrected  by  replacing  t by  r in  (28) 

and  (30).  The  resulting  expression  is  correct  to  first  order;  that  is,  its 

2 

leading-order  error  is  0(e  ).  Since  the  computed  solution  can  be  closer 
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to  (or  further  away  from)  the  correct  nonlinear  solution  than  the 
approximate  analytic  solution,  the  RMS  error  is  an  unknown  mix  of 
errors  in  both  the  computed  and  approximate  solutions. 


4.  2 Finite-difference  model  results 
(Table  I,  cases  1-9;  Figures  1-6) 


An  exact  solution,  to  the  discrete  finite -difference  Eqs.  (5a) 
and  (5b)  can  be  found  by  assuming 

fd(x»y,t)  = sin(Xx)  sin(|i.y)  cos  (ax  + at)  (31) 


where  a and  a are,  in  general,  functions  of  the  nondimens ional  space 
and  time  increments  h and  At.  (We  therefore  assume  that  the  discrete 
and  analytic  results  differ  only  in  the  wavelength  and  phase  speed  of  the 
travelling  component  of  the  box  mode. ) Substituting  in  the  trial  solution 
(31),  we  find  that 

sin(oAt)  = 1/2  hAt 


and 


cos(ah) 


cos  h(\h) 
2-cos(ph) 


For  h « 77  and  At  « 17, 

2 

° = 1 - 12  (2  + ^2)h2  + 0(h4) 


and 

a = 1/2  + 1/48  At2  - l/24(4-2p2-jx4)h2  + 0(At4+h4)  . 

That  this  is  indeed  the  correct  computational  result  has  been  verified  by 
direct  numerical  integration  of  the  finite -difference  equations  (5a)  and 
(5b).  The  resulting  computational  solution  differs  from  (31)  only  due  to 
machine  truncation  error--that  is,  RMS(£’)  = RMS(£c-Cd)  = 0(10-11). 

In  comparison  to  the  analytic  solution  (for  which  a = 1 and  a = 1/2), 
the  wavelength  and  period  are  correct  only  to  second  order  in  space  and 
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time.  For  sufficiently  small  phase  errors  0 (=  a - a,),  it  is  easy  to 

2L  Q 

show  that 

RMS(f')  = [2(1  - costft)]1/2  . 

The  initial  error  growth  rate 

RMS(i|r')  - 0t  + O(0t)3  , (32) 

is  therefore  linear  in  time  with  a slope  given  by 

0 = (l/48)At2  - 1/24(4  - 2p2  - p4)h2 

^ (tr2/3>n"2  - tt2/24(4  - 2p2  - p4)v"2  . (33) 

An  example  of  this  behavior  is  given  in  Figure  (3a,  b)  for  (v,r|,m,  n)  = 
(16\£,  128,  1,  1). 

Note  that  for  the  discrete  finite -difference  solution,  the  spatial  and 
temporal  contributions  to  the  phase  error  0,  being  of  opposite  sign,  tend 
to  offset  each  other.  Because  of  this  compensation  effect,  if  an  optimal 
choice  of  At  and  h is  made,  the  total  error  of  the  finite -difference 
scheme  can  be  made  quite  small  although  the  contributions  to  0 from 
spatial  and  temporal  error  are  individually  large.  This  property  explains 
the  computed  results  as  a function  of  v(“*  v/h)  and  T](=  4ir/At)  in  which 
increasing  q (holding  v fixed),  and  vice- versa,  can  increase,  rather 
than  decrease  the  error  of  the  computed  finite -difference  solution. 
Compare,  for  instance,  Table  1,  cases  2 and  3. 

For  e > 0,  the  values  of  £ on  the  boundary  enter  the  problem 
through  the  nonlinear  terms.  Three  ways  of  fixing  £ ^ have  been 
examined  in  the  context  of  the  finite- difference  model.  They  are  the 
specification  of  the  analytic  value  of  the  vorticity  (£  = £ ),  and  the 

Kreiss  and  Sundstrom/Davies  conditions --see  Section  3.  For  the  non- 
linear box  mode  problems  studied,  the  following  behavior  was  noted. 

4.  2.  1 e = 0. 2 

For  low  to  moderate  Ross  by  number  the  FD  model  is  always  well 
behaved  out  to  at  least  t = 5 periods  independent  of  computational 
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boundary  condition.  When  CE  = Ca  is  the  required  condition,  however, 
there  is  a buildup  of  small-scale  perturbation  vorticity  on  the  western 
wall- -Figure  6a.  (Since  similar  effects  are  noted  in  all  the  models, 
this  buildup  is  presumably  a manifestation  of  the  physical  response  of 
the  system  to  the  presence  of  small-scale  numerical  truncation  error.  ) 
This  accumulation  of  Q'  is  less  rapid  when  the  Kreiss  boundary  condition 
is  applied  everywhere  (Figure  6b)  and  is  nearly  eliminated  when  the 
Sundstrom/Davies  condition  is  invoked  (Figure  6c).  The  RMS  error 
measures  are,  however,  comparably  large--several  tens  of  percent 
after  five  periods --for  all  three  boundary  conditions  (Table  1,  case  6). 

4.  2.  2 e = 0.4 

For  higher  Rossby  number,  the  error  accumulation  to  the  west  is 
much  more  rapid  and  becomes  noticeable  in  even  the  Sundstrom/Davies 
experiments.  In  contrast  to  the  finite -element  and  spectral  models,  how- 
ever, the  finite -difference  scheme  does  not  suffer  catastrophic  numerical 
instability  when  grid- scale  vorticity  begins  to  accumulate.  This  lower 
sensitivity  to  the  presence  of  small-scale  vorticity  is  perhaps  due  to  a 
small  amount  of  (numerical)  dissipation  implicit  in  the  finite -difference 
formalism. 

In  general,  the  RMS  error  quantities  have  a linearly  increasing 
trend  similar  to,  but  somewhat  greater  than,  that  noted  for  e = 0.  Some 
our  approximate  reference  solution. 

4.  3 Finite- element  model  results 
(Table  I,  cases  1-9)  Figures  1-5) 

The  following  functional  dependence  on  the  parameters  v and  q 
has  been  noted  in  the  error  analysis  of  the  finite -element  model  results. 
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r 


4.  3.  1 r)  « 4 v 


For  relatively  coarse  temporal  resolution,  the  RMS  error 

quantities  are  well  described,  as  in  the  finite- difference  and  spectral 

models,  by  the  relation  0 “*  (ff2/3)ri~2  and  consequently  by  an  initial 

_2 

linear  increase  with  time  proportional  to  r]  . This  reflects  the  fact 
that  for  sufficiently  large  v the  error  is  attributable  to  temporal 
truncation  effects  (identical  in  all  three  models). 


4.3.2  r|  =*■  4 \fl  v 

As  r|  increases,  RMS(£')  for  the  finite-element  is  characterized 
by  a large  initial  value  followed  by  a very  slow  linear  increase  there- 
after (Figure  3c).  (RMS(\|r ')- -whose  values  are  perhaps  dictated  by  a 
different  mechanism  than  those  of  RMS(£ ')- -maintains  a linear  trend.  ) 
Apparently,  the  spatial  contribution  to  the  phase  error  0 is  of  a form 

quite  different  than  that  implied  by  (33).  For  one  thing,  we  know  the 

4 • • • 

phase  error  to  be  proportional  to  h for  the  linearized  advective 
equation  [22].  In  parameter  ranges  where  this  diminished  error  growth 
rate  prevails,  the  finite -element  model  may  offer  some  advantages  for 
long-term  calculations;  however,  we  have  not  verified  this  possibility. 

For  e = 0,  the  maximum  pointwise  errors  in  vorticity  tend  to  be 
near  the  western  boundary,  but  there  is  very  little  preferential 
accumulation  of  small-scale  vorticity  there.  With  nonlinearity  (e  >0), 
the  situation  is  qualitatively  different  in  the  following  way(s). 

4.  3. 3 e = 0. 2 

With  € = 0.2,  an  initial  eastern  boundary  layer  is  generally 
observed  in  the  field  of  Q'.  This  boundary  layer  eventually  disappears, 
to  be  replaced  by  an  accumulation  of  perturbation  vorticity  on  the 
western  wall,  as  in  the  FD  and  PS  simulations  (Figure  5c).  In  other 


L 
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respects,  the  solutions  bear  some  resemblance  to  those  for  e = 0. 

RMS(C')  again  shows  evidence- -for  certain  values  of  v and  t)--of 
levelling  off  with  time  after  an  initially  large  increase,  and  i(r ' some- 
times resembles  a box  mode  (out  of  phase  with  the  reference  solution). 

4.  3.4  e = 0.4 

With  stronger  nonlinearity,  perturbation  vorticity  on  the  grid-point 
scale  collects  first  on  the  western  wall  and  then  in  the  center  of  the 
domain  (perhaps  as  a numerical  response  to  insufficient  resolution  of 
the  narrow  wall  layers  of  £').  Once  this  stage  is  reached,  the  solution 
becomes  numerically  unstable,  typically  after  about  5 periods  (Table  I, 
case  7). 

This  catastrophic  effect  of  small-scale  vorticity  accumulation  is 
reached  in  less  than  a period  for  e = 0.  8. 

4.  4 Pseudospectral  model  results 
(Table  I,  case  1-9;  Figures  1-5,7) 

For  linear  box  modes,  for  which  we  have  the  analytic  solution,  the 
spectral  model  shows  three  distinct  types  of  behavior  corresponding  to 
different  regimes  in  the  space  of  the  nondimensional  computational 
parameters. 

4.4.1  v 4^2,  r)  < 128  (t)  « 1 6 /2  v ) 

Quite  a large  range  of  v exists  for  which  spatial  truncation  errors 

are  totally  insignificant  in  comparison  to  those  arising  from  time- 

differencing.  For  this  range  of  parameters,  the  RMS  quantities  grow 

linearly  in  time  (Figures  3e,f)  and  can  be  quantitatively  explained  by  the 

simple  phase  error  analysis  of  Section  4.  2 if,  in  addition,  the  assumption 

is  made  that  h “■  0.  There  is  no  apparent  buildup  of  perturbation 

vorticity  at  scales  other  than  those  of  the  box  modes  themselves.  Since 

the  computational  errors  are  due  to  time -differencing  alone,  they  are 
-2 

proportional  to  q --see  Table  I,  cases  1 and  2. 
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4.4.2  -^svs4^,  n < 128  (ri  - 16^v) 

For  v “*  4 and  standard  values  of  r|(=  64  or  128),  spatial  error 
becomes  noticeable.  The  growth  of  the  RMS  quantities  is  no  longer 
strictly  linear.  In  fact,  the  temporal  variation  of  RMS(£')  begins  to 
assume  characteristics  noted  (over  a larger  range  of  v and  r|)  in  the 
finite -element  model:  a large  initial  error  growth,  followed  by  a 
relatively  slow  increase  with  time.  Indeed,  it  is  interesting  to  note 
that  the  spectral  model  can  show  the  effects  of  compensating  space  and 
time -differencing  errors  which  is  a general  property  of  finite-difference 
model.  Something  of  this  kind  is  clearly  happening  in  the  spectral  model 
when,  for  instance,  RMS(|')  decreases  when  the  box  mode  numbers 
m = n are  increased  from  two  to  three  at  constant  v and  r]  (Table  I, 
cases  4 and  5).  For  these  values  of  v,  perturbation  vorticity  does  appear 
to  collect  on  the  western  wall,  perhaps  in  scales  much  smaller  than  those 
of  the  original  box  modes. 

4.4.3  v < Z'JI,  rj  * 128  (rj»  16^2  v) 

As  expected,  for  extremely  small  values  of  the  resolution  para- 
meter v,  substantial  spatial  error  results.  The  RMS  error  quantities 
once  again  grow  in  a quasilinear  fashion.  The  perturbation  vorticity 
field  is  now  dominated  by  a narrow  layer  on  the  western  wall.  The 
amplitude  of  this  feature  is  sufficiently  large  (after  10  periods)  as  to 
contribute  recognizably  to  the  total  vorticity  field. 

For  the  nonlinear  box  modes,  behavior  of  the  computational  system 
depends  sensitively  on  e in  the  following  manner. 

4.4.4  e = 0.  2 

For  e <0.2,  the  spectral  model  is  well  behaved  out  to  t — 5 

2 

periods.  By  this  time,  however,  integrated  |v£|  has  begun  to  increase 
rapidly.  Significantly  longer  integrations  could  therefore  be  expected 
to  suffer  eventual  computational  instability.  Even  at  this  level  of 
nonlinearity,  increasing  v and/or  t)  (over  the  range  tested:  v s 8 -Ji, 

T)  £ 128)  does  nothing  to  improve  the  error  measures  (Table  I,  cases  6, 
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8,9).  The  manifestation  of  error  growth  is  a very  definite  preferential 
accumulation  of  perturbation  vorticity  in  narrow  layers  adjacent  to  the 
western  wall  of  the  domain  (Figure  5e).  This  accumulation  of 
perturbation  vorticity  may  ultimately  result  from  local  numerical 
truncation  errors  which  are  propagated  to  the  west  where,  in  the 
absence  of  dissipation,  they  collect  in  a narrow  boundary  layer. 

4.4.5  e = 0.  4 

The  results  for  e = 0.  4 are  much  more  catastrophic,  with  per- 
turbation vorticity  collecting  so  quickly  on  the  western  wall  that  locally 
intense  gradients  of  vorticity  grow  to  destroy  the  calculation  after  only 
2.  5 periods  (Table  I,  case  7).  This  behavior  is  once  again  independent 
of  v and  r\.  When  the  calculations  go  bad,  they  do  so  very  quickly;  pre- 
sumably the  computed  fields  are  still  quite  accurate  up  to  the  instant  of 
catastrophic  failure.  (This  essentially  instantaneous  instability  is  a 
feature  of  the  FE  model  also.  ) In  a related  calculation,  it  has  been  shown 
that  the  useful  integration  of  the  spectral  model  can  be  prolonged  to  t = 5 
periods  (and  beyond)  by  periodically  filtering  the  vorticity  field  by  setting 

bk  (filtered)  = f f 
nm  n m nm 

(see  Section  3.  3)  where  the  spectral  filter 

f = 1.0  - exp{-J(f(N^  - n^)} 

and  is  adjusted  so  that  £ is  smoothed  only  at  the  highest  wave- 
numbers.  By  comparing  the  filtered  and  unfiltered  results,  it  is  known 
that  such  filtering  does  not  affect  the  large-scale  features  of  the 
circulation  and  that  the  two  streamfunction  fields  (up  to  the  moment  of 
instability  in  the  unfiltered  calculation)  are  virtually  identical  (Figure  7). 

Lastly,  it  is  important  to  note  that  the  RMS  error  quantities  for 
the  nonlinear  box  mode  problems  are  nearly  independent  of  v and  T]. 

If  there  is  nothing  idiosyncratic  about  these  problems,  then  we  must 
conclude  that  the  largest  contribution  to  the  RMS  error  fields  comes 
from  the  uncertainty  in  the  exact  analytic  solution  to  the  nonlinear  box 
mode  problem. 
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4.  5 Intercomparison 

Results  from  the  linear  box  mode  tests  (see  Table  I)  demonstrate 
that  the  finite -element  (with  N = 33)  and  pseudospectral  (with  N = 17) 
models  are  comparably  accurate  over  the  range  of  v and  17  studied. 

(The  spectral  model  is,  however,  somewhat  more  efficient-- Table  II.  ) 

And,  even  though  the  finite -difference  is  by  far  the  least  accurate  model, 
a phase  error  analysis  of  the  FD  model  results  shows  that  errors  can  be 
minimized  for  optimal  choices  of  v and  r).  Since  these  optimal  para- 
meters are  functions  of  the  time  and  space  scales  of  the  problem,  however, 
this  property  will  be  of  questionable  value  in  more  general  problems 
characterized  by  multiple  time  and  space  scales.  Even  if  a degree  of 
compensation  could  be  guaranteed  in  a specific  problem,  errors  get 
smaller  only  if  v and  77  are  increased  in  the  same  ratio.  A fourfold 
decrease  in  the  RMS  errors  would  therefore  require  17  and  v to  be 
simultaneously  increased  by  a factor  of  two,  with  a resulting  increase  in 
computational  work  of  a factor  of  eight.  Jn  contrast,  the  spectral  model 
(and  to  a lesser  extent  the  FE  model)  generally  require  only  17  to  be  in- 
creased--say  by  two,  for  a four-fold  reduction  in  error- -because  of  their 
much  greater  spatial  accuracy. 

In  the  case  of  the  nonlinear  box  modes,  interpretation  of  the  results 
is  complicated  by  the  fact  that  we  have  only  a perturbation  solution  with 
which  to  compare  the  computational  results.  Consequently,  our  error 
measures --such  as  RMS(\|r'),  etc. --reflect  three  sources  of  error:  spatial 
and  temporal  truncation  errors,  and  the  error  associated  with  not  knowing 
the  exact  analytic  solution.  The  RMS  quantities  listed  in  Table  I for  the 
nonlinear  box  mode  tests  cannot  be  used  as  direct  measures  of  model 
performance. 

The  essential  qualitative  distinction  that  can  be  made  between  the 
results  of  the  three  models  for  e >0  is  that  the  finite -difference  model, 
though  presumably  less  accurate,  appears  not  to  be  susceptible  to 
catastrophic  numerical  instability  when  small-scale  error  fields  are 
present.  Under  these  conditions,  the  FO  spatial  truncation  error  is,  how- 
ever, formally  quite  large.  A nonlinear  FO  solution  will  therefore  become 
invalid  after  only  a short  period  of  time  even  though  a stable  calculation 
can  be  maintained  for  a much  longer  time. 
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5.  FORCED  NONLINEAR  BOX  MODE  TESTS 
5.  1 Formulation 

One  means  of  avoiding  the  complications  associated  with  having 
only  a perturbation  solution  to  the  nonlinear  box  mode  problem  is  to 
consider  the  analogous  forced  problem,  that  is  to  seek  solutions  to  the 
inhomogeneous  equation 

~ VZf  + J(t,  VZf)  + \|tx  = F(x,  y,t)  (3 

where  F is  some  suitably  chosen  forcing  function.  As  before,  t|r  is 
required  to  vanish  on  2 . In  particular,  we  wish  to  examine  solutions 
with  spatial  and  temporal  characteristics  similar  to  those  of  the  linear 
box  modes.  Accordingly,  set 


i|r(x,  y,  t)  = sinx  siny  cos(ax  + by  + ct) 


0 £ x,  y s n 


where  a,  b,  and  c are  arbitrary  constants  which  determine  the  wave- 

2 2 1 /2 

length,  period  and  phase  speed  of  the  forced  mode  (27r/(a  +b  ) ' , 2tt/c, 

2 2 1 / 2 

and  c/(a  +b  ) ' , respectively).  This  will  be  a solution  to  (34)  so  longas 


F(x,y,t)  = j-^  V2  + J(i|r,eV2)  + -g- J 


which  will  in  general  be  nonvanishing,  as  will  J(i|r,  ev  y).  Given  specific 
values  of  parameters  a,  b,  c and  e,  the  functional  form  of  F can 
therefore  be  directly  calculated.  For  the  following  tests,  the  Rossby 
number  has  been  fixed  at  e = 0.2.  An  examination  of  higher  e behavior 
is  reserved  for  the  open  boundary  calculations  of  Section  7. 

5.  2 Finite -difference  model  results 
(Table  I,  cases  10-12;  Figures  8-10) 

Table  I shows  the  RMS  error  measures  for  the  finite -difference 
model  after  two  periods  for  a variety  of  values  of  v and  t)  . The  results 
indicate  that  the  FD  error  norms  are  in  general  somewhat  smaller  than 
those  for  the  linear  box  mode  problems  with  comparable  nondimens ional 
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parameters.  Compare  Figures  (3a,  b)  and  (10a,  b),  for  instance.  In 
addition,  a partial  compensation  between  spatial  and  temporal  errors 
once  again  exists  so  that  the  RMS  errors  (as  in  the  linear  box  mode  cases) 
need  not  decrease  with  increasing  rj  and  v (Table  I,  cases  11  and  12). 

As  in  the  nonlinear  box  mode  problems,  perturbation  vorticity  tends  to 
collect  on  the  western  boundary  (Figure  9a).  This  appears  to  be  a quite 
general  property  of  all  the  simulations  when  e >0,  no  matter  what  the 
orientation  of  the  forced  mode. 


5.  3 Finite -element  model  results 
(Table  I,  cases  10-12;  Figures  8-10) 


The  FE  model  behaves  similarly,  yielding  very  accurate  and  stable 
solutions  for  a range  of  parameters  (Table  I,  cases  10-12).  For  instance, 
with  (v,rj,a,b)  = (32,128,1/'^,  l/*j2)  the  RMS  errors  are  0 (1-25?)  after 
two  periods.  As  in  the  FD,  and  as  we  shall  see,  the  spectral  computations, 
the  RMS  errors  given  by  the  FE  solution  to  the  nonlinear  forced  box  mode 
problem  are  typically  no  less,  and  very  often  several  times  smaller,  than 
the  errors  noted  for  the  linear  unforced  box  mode  tests  with  comparable 
resolution.  Figures  (3c,  d)  and  (10c,  d)  give  an  example  of  this  behavior. 
(Note,  also  that  the  character  of  the  RMS  error  curves  seems  to  be  modi- 
fied by  the  forcing  such  that  RMS(£')  is  a quasilinear  function  of  time  over 
the  range  of  parameters  examined  here.  ) The  perturbation  fields  associated 
with  the  forced  problems,  although  small  in  amplitude,  are  still  character- 
ized by  small-scale,  westward-trapped  Q'  and  large-scale  f ' patterns 
(Figure  9 )• 


5.4  Pseudospectral  model  results 
(Table  I,  cases  10-12;  Figures  8-10) 

These  calculations  were  all  performed  with  very  high  spatial 
resolution  (v  £ 16);  consequently,  insignificant  spatial  truncation  error 
is  expected.  With  a diagonally  propagating  mode  (a  = b = 1 /\[Z)  and 
r)  = 128,  the  RMS  errors  are  in  fact  very  small,  being  no  more  than 
0.  5 percent  after  5 periods  of  integration.  Reference  to  Table  I and 
Figures  10e,f  demonstrates,  however,  that  not  only  are  the  RMS  errors 


no  longer  strictly  proportional  to  r|”  --as  they  were  in  the  unforced 
case  with  sufficient  spatial  resolution- -but  the  error  trends  are  not 
linear  in  time.  Although  RMS(£')  increaees  monotonically,  RMS(f') 
seems  to  vary  quasiperiodically  with  little  superimposed  trend 
(Figures  10e,f). 

As  in  the  FD  and  FE  simulations,  the  PS  results  for  the  forced 
nonlinear  box  mode  problems  show  a similar  tendency  for  small-scale 
Q'  to  accumulate  at  the  western  edge  of  the  domain  (Figure  9e).  The 
rate  of  this  error  accumulation  might  plausibly  be  thought  to  increase 
dramatically  with  e,  as  in  the  nonlinear  box  mode  problems;  however, 
this  hypothesis  was  not  tested. 

None  of  these  conclusions  depend  sensitively  on  the  direction  of 
propagation  of  the  forced  mode. 

5.  5 Inter  comparison 

All  three  models  are  capable  of  delivering  an  accurate  and  stable 
solution  to  the  forced  nonlinear  box  mode  problem  for  the  computational 
parameters  considered  here.  In  all  cases,  the  observed  RMS  errors 
are  less  than  or  equal  to  those  noted  in  the  linear  box  mode  cases  (for 
comparable  v and  r)).  It  is  quite  likely  that  this  reduction  in  numerical 
error,  despite  going  to  a problem  with  nontrivial  nonlinearities,  is  in 
some  sense  associated  with  a "locking  in"  of  the  numerical  solution  to 
the  applied  forcing.  The  quasioscillatory  nature  of  ^he  resultant  RMS(|') 
error  curves--see  especially  Figure  10f-- argues  for  such  3.  process. 

We  will  return  to  this  possibility  in  our  discussion  of  future  work. 


Consider  solving  Eq.  (3)  with  e = 0 in  a domain  characterized  by- 
three  closed  boundaries  (on  the  east,  south  and  west),  but  which  is  open 
on  its  northern  edge.  If,  in  particular,  we  let  the  boundary  conditions  on 
2 be 


x=0,x„;  y=  0 


sin  y sin 


in(^) 


cos(x  + t/2) 


then  the  solution  is 


ijf(x,y,t)  = sin  sin|-^|  cos (x  + t/2) 


22  21/2 

provided  that  x^  =yB=(mir  + y ) ' • A special  case  of  these  north 
wall  forced  modes  are  the  linear  box  modes,  Eq.  (28),  which  result  when 
y = n7T,  n an  integer,  so  that  the  streamfunction  vanishes  on  all  four  walls. 

Under  appropriate  forcing  conditions,  trapped  solutions  also  exist 
to  the  north  wall  forced  problem.  One  simple  solution  is 


(r(x,  y,t)  = sin  sinh  ) cos(x  + t/2) 


where 


= yR  = (m  jt  - Y ) 


2.1/2 


The  appropriate  boundary  conditions  for  this  problem  are  clearly  the 
value  of  if  evaluated  along  the  boundary  2.  Since  these  solutions  decay 
away  from  the  north  wall,  (39)  is  referred  to  as  a trapped  case,  and  (38) 
as  a propagating  case.  Note  that  for  propagation  (trapping)  xfi  > mir 
(x_  < mir).  A more  general  description  of  the  north  wall  forcing  problem 
and  its  relation  to  meander-induced  forcing  has  been  given  by  Harrison 
[31]. 
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Two  particular  north  wall  forcing  problems  have  been  considered, 
one  each  of  the  propagating  and  trapped  varieties.  The  former,  which 
we  will  refer  to  as  case  1,  corresponds  to  Eq.  (38)  with  m = 1 and 
■y  = 9.  36.  The  latter  experiment,  or  case  2,  is  the  trapped  solution  (39) 
with  m - 2 and  y = 4.0. 

In  addition  to  these  linear  solutions,  the  corresponding  nonlinear 
case  1 solution- -that  is,  the  solution  to  Eqs.  (3)  and  (37)  with  m = 1, 
y = 9.  36  and  e > 0--has  been  sought.  With  e = 0,  problems  1 and  2 
are  similar  in  practice  to  the  linear  box  mode  problems,  except  that 
assumes  some  time -dependent  distribution  of  values  along  the  northern 
boundary.  With  nonlinearity  present,  however,  a set  of  vorticity 
boundary  conditions  must  also  be  specified  along  y = y^.  In  the  FE  and 
PS  models,  vorticity  is  given  its  analytic  value  on  boundary  points 
characterized  by  inflow.  In  addition  to  this,  the  FD  model  constrains 
the  values  of  vorticity  on  outflow  by  the  Sundstrom/Davies  condition- -see 
Section  3.  1.  For  the  purposes  of  this  study,  solid  boundaries  will  be 
treated  as  points  of  outflow.  Since  an  exact  solution  to  the  resulting  non- 
linear problem  is  not  available,  we  use  (38)  as  our  reference  solution 
even  when  e > 0.  This  leads  to  a problem  in  interpreting  the  RMS  error 
quantities  similar  to  that  encountered  in  the  nonlinear  box  mode  tests 
where  only  a perturbation  solution  was  available. 

6.  2 Finite -difference  model  results 
(Table  I,  cases  13-21;  Figures  11-19) 

Not  only  are  the  methods  of  solution  nearly  identical,  as  noted 
above,  but  the  results  of  the  cases  1 and  2 linear  north  wall  forced  mode 
problems  are  themselves  quite  similar  to  those  of  the  linear  box  mode 
tests  (Section  4).  The  following  v and  ri  dependencies  were  noted  in 
the  propagating  and  trapped  cases,  respectively. 

6.  2.  1 Case  1 (propagating) 


The  case  1 north  wall  forced  mode  closely  resembles  the  linear 
box  mode  with  m = n = 3.  The  results  for  these  two  FD  test  problems 
are  comparable  with  few  exceptions.  The  error  fields  are  once  again 
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attributable  to  a simple  phase  difference  between  the  discrete  and 
analytic  solutions  (Figure  11a,  b).  RMS(ji')  and  RMS(£')  grow  linearly 
in  time  over  the  five  period  integrations  (Figure  13a,  b).  Because  of  the 
highly  structured  nature  of  these  solutions  and  the  second-order  spatial 
accuracy  of  the  model,  the  RMS  FD  errors  tend  to  be  quite  large  in  both 
although  their  variations  with  v and  rj  (Table  I)  suggest  that  it  is 
possible  to  reduce  the  RMS  error  measures  for  suitable  choices  of  the 
computational  parameters.  One  striking  dissimilarity  between  the  box 
mode  and  north  wall  forced  mode  results  is  that,  of  the  three  error 
measures  tabulated,  only  the  normalized  error  in  integrated  energy  is 
greatest  for  the  north  wall  problem.  This  presumably  reflects  the 
modified  nature  of  the  north  wall  forced  problem  in  which  energy  is 
exchanged  between  the  interior  solution  and  exterior  environment  at  a 
rate  determined  by  the  imposed  values  of  \|r^  and  . The  change  in 
energy  level  over  the  course  of  the  5 period  simulation  is,  however, 
always  fractionally  quite  small,  being  no  more  than  0 (10$)  even  for  the 
poorest  resolved  experiment. 

6.  2.  2 Case  2 (trapped) 

Although  the  trapping  scale  y * in  equation  (39)  can  be  quite  small, 
it  has  been  chosen  in  case  2 to  give  only  moderate  structure  to  the 
solution.  Consequently,  the  case  2 solution  is  actually  somewhat  less 
structured  than  the  propagating  case  and,  as  such,  the  trapped  FD  results 
are  more  accurate  by  roughly  a factor  of  four  (Table  I,  cases  14  and  17, 
for  instance).  The  curves  of  RMS(£')  tend  to  rise  to  an  initial  moderate 
value,  but  grow  much  more  slowly  thereafter  (Figure  16a)--a  rattern 
noted  over  varying  ranges  of  the  parameters  in  the  PS  and  FE  linear  box 
mode  experiments.  (RMS(^')  tends  to  remain  quasilinear  in  time  though 
it  also  has  large  excursions  about  the  linear  trend.  ) For  the  case  2 
simulations,  the  perturbation  streamfunction  and  vorticity  fields  can  no 
longer  be  explained  by  a simple  phase  difference  between  computed  and 
analytic  fields.  Instead,  perturbation  vorticity  tends  to  accumulate  in 
the  northeast  and  northwest  corners  of  the  domain  (Figure  15a). 

For  e = 0.2,  the  character  of  the  FD  solution  is  changed  dramatic- 
ally. Early  in  the  simulation,  the  fields  of  perturbation  streamfunction 


and  vorticity  are  characterized  by  a 6x2  celled  structure.  In  fact, 
this  feature  appears  in  the  FE  and  PS  simulations  as  well  (Figure  18). 

Since  in  this  problem  our  reference  solution  is  the  linear  forced  mode 
(38),  this  structured  perturbation  vorticity  field  must  partially  reflect 
the  nonlinear  correction  to  (38)  for  e = 0.  2.  At  later  stages  in  the 
simulation,  the  perturbation  vorticity  generated  during  the  calculation 
winds  up  on  the  western  boundary  where  it  forms  six  narrow  cells  of 
vorticity  adjacent  to  the  walls  (Figure  19).  For  the  FD  model,  these 
layers  of  vorticity  dominate  the  field  of  total  vorticity  by  the  end  of  two 
periods  of  integration.  Although  this  graininess  in  the  FD  results  does 
not  typically  lead  to  numerical  instability,  such  Q fields  are  not  accurately 
resolvable  on  the  FD  grid.  For  e = 0.4,  this  state  of  an  unresolvable 
vorticity  field  is  reached  after  one  period. 

As  in  the  nonlinear  box  mode  problem,  the  RMS  errors  are  not 
calculated  with  respect  to  the  exact  analytic  solution. 

6.  3 Finite -element  model  results 
(Table  I,  cases  13-21;  Figures  11-19) 

As  noted  for  the  FD  model,  the  qualitative  characteristics  of  the 
case  1 and  case  2 FE  computational  results  are  very  different, 
particularly  in  the  form  of  the  i|r 1 and  £'  fields  (Figures  12c,  d and  15c,  d) 
and  the  associated  RMS  error  curves  (Figures  13c,  d and  16c,  d). 

6.  3.  1 Case  1 (propagating) 

As  expected  from  their  structural  similarity,  the  3x3  linear 
box  mode  and  case  1 north  wall  forced  mode  problems  both  give  RMS 
errors  of  0 (5$)  after  five  periods  for  (v,r|)  = (0(10),  64)--see  Table  I, 
cases  5 and  14.  (This  is  to  be  compared  with  errors  of  30%  and  greater 
for  the  comparable  FD  problems.  ) The  perturbation  fields  are  them- 
selves box  mode -like,  indicating  a simple  phase  error  relationship 
between  computed  and  analytic  solutions.  Depending  on  v and  r|,  the 
RMS(C')  error  curves  may  have  either  of  the  two  forms  noted  in  the 
linear  box  mode  calculations.  For  the  pivotal  resolution  (v,T})=  32/3,64), 
both  RMS(if')  and  RMS(£')  increase  quasilinearly  (Figure  13c,  d). 
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While  RMS(C')  decreases  with  increasing  v and  q,  RMS(f')  does  not, 
perhaps  indicating  some  degree  of  compensation  between  the  contributions 
to  ijr  of  spatial  and  temporal  errors. 

6.  3.  2 Case  2 (trapped) 

For  case  2,  the  fields  of  perturbation  streamfunction  and  vorticity 
are  not  simply  related  to  the  analytic  form  of  the  trapped  modes.  While 
is  dominated  by  lateral  scales  larger  than  those  making  up  the  trapped 
modes,  Q'  appears  at  highest  wave  numbers  and  is  eventually  trapped  in 
layers  on  both  the  eastern  and  western  walls  (Figure  15c,  d).  The  RMS 
errors  increase  only  slowly  with  time  (Figure  16c,  d)  and  decrease 
uniformly  with  increasing  v and  q around  the  point  of  pivotal  resolution 
(Table  I).  Typical  error  after  five  periods  are  0 (2-3$)  for  RMS(f') 
and  RMS(G')  at  (v,q ) = (64/3, 64). 

Stable,  but  noisy,  solutions  can  be  obtained  out  to  five  periods  for 
the  case  1 problem  with  e = 0.  2.  The  final  fields  for  the  pivotal  case-- 
(v , T))  = (32/3,  64)--are  given  in  Figures  17c,  d and  19c,  d.  Note  the 
preferential  accumulation  of  grid-point  scale  vorticity  in  the  western- 
most basin  whose  contribution  to  the  total  field  Q is  evident  as  early  as 
t = 2 periods.  Although  gridpoint  variations  in  the  west  undoubtedly 
contribute  to  sizeable  numerical  discretization  errors  in  this  region,  it 
must  be  added  that  the  resulting  total  streamfunction  field  i|f  appears  to 
be  only  locally  affected,  retaining  over  most  of  the  basin  the  expected 
mode-like  form.  At  6=  0.4,  the  small-scale  vorticity  accumulates 
much  more  catastrophically,  leading  to  numerical  instabilities  for  ts  5 
periods  (Table  I,  case  21). 

6.4  Pseudospectral  model  results 
(Table  I,  cases  13-21;  Figures  11-20) 

The  results  of  the  Chebyshev  spectral  model  are  both  qualitatively 
and  quantitatively  similar  to  those  of  the  finite -element  code.  This  is 
perhaps  due  to  their  higher  order  spatial  accuracy. 
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6.4.  1 Case  1 (propagating) 

For  the  values  of  v used  in  the  case  1 tests,  the  spectral  model 

has  negligible  spatial  differencing  error.  The  perturbation  fields  are 

_ 2 

therefore  attributable  to  a phase  error  proportional  to  r)  and  not  at 
all  to  v (Figures  lie,  f and  12e,f;  Table  I).  For  (v,ri)  = (32/3,64),  the 
resulting  errors  are  0 (5$)- -comparable  to  (much  less  than)  the  FE  (FD) 
error  norms.  RMS(i)1)  and  RMS(C')  increase  linearly  in  time 
(Figure  13e,f). 

6.4.2  Case  2 (trapped) 

Although  i|t ' and  Q 1 are  no  longer  characterized  by  a simple  phase 
relationship  to  \|r  and  Q , but  by  basin-scale  and  grid-scale  features, 

d a 7 

respectively  (Figures  15e,f),  their  RMS  values  at  t = 5 periods  are  again 

_ 2 

approximately  proportional  to  T)  (Table  I,  cases  17  and  18).  Fora 
pivotal  resolution  of  (v,r))=( 64/3,64),  RMS(i|t ')  “"RMStC ')  = 0(1$).  After 
an  initial  period  of  rapid  error  growth,  the  RMS  norms  are  oscillatory  in 
time  with  a much  slower  linear  trend. 

When  we  set  e = 0.  2 and  seek  nonlinear  solutions  to  the  case  1 north 
wall  forced  mode  problem,  the  spectral  model  goes  unstable  as  early  as 
t = 2 periods.  This  behavior  coincides  with  that  of  the  FE  model,  but  at 
slightly  higher  e.  The  origin  of  this  instability  is  the  ultrathin  western 
boundary  layer  of  vorticity  which  develops  in  all  the  models  by  this  time 
as  a result  of  accumulated  computational  error  (Figure  19e).  The 
Chebyshev  spectral  model  is  especially  sensitive  to  such  a feature  of  the 
solution  because  of  its  nonuniform  distribution  of  resolution  which  favors 
the  domain  edges  over  the  interior.  For  reasonably  narrow  layers,  the 
Chebyshev  model  will  give  a much  more  accurate  representation  than 
either  the  FD  or  FE  models.  However,  the  accuracy  of  the  spectral 
model  works  to  its  own  disadvantage  for  layers  which  are  much  thinner 
than  the  collocation  grid  spacing.  This  relationship  between  the  accuracy 
and  stability  of  FD  and  PS  approximation  has  also  been  noted  in 
integrations  of  one -dimensional  viscous  transport  problems  [15]. 

The  most  straightforward  way  to  avoid  eventual  computational 
instability  is  to  suppress  the  generation  of  perturbation  vorticity  and  its 
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accumulation  on  the  western  boundary.  This  could  be  done  by  greatly 
increasing  the  spatial  and  temporal  resolution  of  the  spectral  model. 

This,  however,  is  inefficient.  Quite  good  results  have  been  attained  by 
periodic  filtering  of  the  vorticity  field  to  preferentially  remove  numeric- 
ally generated  vorticity  errors  which,  it  is  assumed,  are  of  much 
smaller  scale  than  the  solution  we  seek.  Such  an  exponentially  tapered 
spectral  filter--see  Section  4.4--has  been  applied  to  the  pivotal  non- 
linear north  wall  forced  experiment--(v,q,  c)  = (32/3,  64,  0.  2)--with  the 
results  shown  in  Figure  20.  After  1.  5 periods,  the  streamf unction  field 
(which  reflects  the  large-scale  component  of  the  flow)  of  both  the  unfiltered 
and  filtered  simulations  are  virtually  identical.  On  the  small  scales, 
localized  layers  of  £'  exists;  however,  they  are  of  much  smaller 
amplitude  and  not  as  thin  in  the  filtered  calculation.  Subsequently,  the 
unfiltered  run  goes  unstable  at  t “■  2 periods.  With  the  filter,  the 
calculation  proceeds  to  t = 5 periods  with  a resulting  field  which  is  at 
least  as  smooth  as  the  FD  and  FE  simulations  (Figure  17).  One  potentially 
undesirable  side -effect  of  the  vorticity  filter  is  that  a much  greater 
fraction  of  the  initial  energy  has  been  lost  over  the  course  of  the  in- 
tegration (Table  I,  case  19).  A more  scale -selective  filter  might  avoid 
this  problem. 


6.  5 Intercomparison 


The  inter  comparative  statements  that  can  be  made  on  the  basis  of 
the  linear  north  wall  forcing  problems  do  not  differ  substantially  from 
those  made  in  connection  with  the  linear  box  modes  - -Section  4.  5.  In 
general,  both  the  FE  and  PS  models,  by  virtue  of  their  superior  spatial 
accuracy,  are  notably  more  accurate  than  the  FD  model.  The  latter, 
however,  can  have  small  overall  discretization  error  when  values  of  v 
and  T)  are  chosen  to  insure  partial  compensation  between  spatial  and 
temporal  errors.  Errors  in  the  PS  simulations  are  strictly  proportional 
to  At^.  This  is  also  approximately  true  for  the  FE  code  in  which,  how- 
ever, there  is  also  a somewhat  complicated  parametric  dependence  on 
the  form  of  the  RMS  error  growth  curves.  For  case  1 (box  mode-like 
solution),  and  £'  themselves  resemble  the  analytic  result  for  all 
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three  models.  In  case  2,  f and  £'  collect  at  low  and  high  wavenumbers 
respectively,  independent  of  numerical  technique. 

When  c > 0,  the  stability  of  the  various  models  seems  to  be  in- 
versely related  to  the  formal  accuracy  of  the  numerical  technique 
involved.  While  the  FD  formalism  is  stable  for  all  the  tests  we  have 
conducted,  the  FE  model  develops  instabilities  for  some  values  of  the 
parameter  set  (v,r|,e).  The  spectral  result  becomes  unbounded  much 
earlier  than  the  FE  calculation,  but  can  be  stabilized  for  integrations  of 
moderate  length  by  spectrally  filtering  the  vorticity  field  to  delay  the 
accumulation  of  vorticity  at  unresolvable  scales. 
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7.  LINEAR  AND  NONLINEAR  ROSSBY  WAVE  TESTS 


7.  1 Formulation 


The  advected  Rossby  wave 

t(x,y,t)  = - yy  + sin(kx  + ty  + out) 


where 

a)  = k(l  - e-y)  (41) 

and 

k2  + |2  = 1 

is  a solution  to  both  the  linear  and  nonlinear  vorticity  equation  (3).  As 
for  the  box  modes,  we  have  scaled  with  respect  to  d = (27t)“  1 times  the 
wavelength  of  the  travelling  wave;  uQ,  a characteristic  particle 
velocity;  and  the  time  scale  (Pd)'1.  In  the  resulting  nondimens ional 
system,  the  wavelength  is  2jt,  and  the  basin  size  x„  = tta  where  A is 

D 

the  number  of  half  wavelengths  per  box  width  (a  measure  of  the  structure 
of  the  solution).  Theoretically  it  is  known  that  Rossby  waves  are 
individually  unstable  to  small  perturbations  [32]  with  an  e-folding  time 
proportional  to  (e)  . This  growth  time  scale  is  comparable  in  all  cases  to 

the  entire  duration  of  the  experiment.  Due  to  the  absence  of  large -amplitude 
perturbations  (or  "noise'')  that  can  efficiently  extract  energy  from  the 
primary  wave,  it  is  unlikely  that  purely  physical  instabilities- -as  opposed 
to  computational  ones --play  a role  in  the  following  results.  The  reader 
should  note  that  the  nonlinearity  of  these  model  problems  is  trivial  (that 
is,  self- cancelling)  when  y = 0. 

7.  2 Finite -difference  model  results 
(Table  I,  cases  22-37;  Figures  21-23,  25-27) 

The  results  for  one  linear  Rossby  wave  experiment  in  which 
(v,T),e)  = (32/3.5,  128,  0)  are  listed  in  Table  I,  case  22.  Variations  in 
the  RMS  error  measures  as  a function  of  v and  t)  did  not  differ  from  the 
comparable  dependencies  noted  for  the  linear  box  mode  and  linear  north 
wall  forced  modes --see  Sections  4.  5 and  6.  5- -and  hence  will  not  be 
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reiterated  here.  The  one  substantive  difference  between  these  and  the 
other  linear  problems  is  that  RMS(f')  and  RMS(C')  are  not  strictly 
linear  in  time  but  appear  to  be  levelling  off  at  t = 5 periods. 

With  no  mean  flow  and  moderate  nonlinearity  (e,  k,f  , y)  = (0.  4, 

3/-/T3,  2/VT3,  0),  the  RMS  errors  of  the  FD  simulations  are  characterized 
by  very  small  temporal,  relative  to  spatial,  errors.  At  the  pivotal 
resolution  (v,T])  = (32/3.5,64),  RMS(f  1 ) “■  RMS(£ ')  •*  NDIF(NRG)  = 0 (14^) 
after  five  periods.  The  manifestation  of  error  in  these  simulations  are 
perturbation  fields  closely  resembling  box  modes  which  begin  to  destroy 
the  plane  wave  nature  of  the  solution  after  a few  periods  (Figures  21a,  b 
and  22a,  b).  This  form  for  the  error  fields  appears  to  be  independent  of 
v and  T].  Their  amplitude,  as  stated  previously,  depends  sensitively  on 
v but  not  on  r)  for  those  values  considered  here.  (In  addition,  other 
simulations --cases  28,  29  and  31- -show  that  RMS(|')  and  RMS(£') 
also  depend  on  the  orientation  of  the  reference  wave  so  that  cancellation 
of  time  and  space  errors  can  sometimes  occur.  ) The  RMS  error 
measures  typically  grow  quasilinearly,  though  in  some  cases  there  is  a 
tendency  for  the  rate  of  error  growth  to  slow  towards  the  end  of  the 
simulation.  Lastly,  the  change  in  the  integrated  kinetic  energy  of  the 
system  is  0 (-IO56),  somewhat  larger  than  that  observed  in  either  the 
nonlinear  box  mode  or  nonlinear  north  wall  forced  mode  problems  for 
comparable  v and  r). 

For  e = 0.8,  but  still  with  (v,ri,k,  I,  y)  = (32/3.  5, 64,  3/VT3,  2/^13,  0), 
the  same  qualitative  remarks  apply.  The  field  of  £'  does,  however, 
begin  to  show  some  noticeable  grid- scale  variability  in  comparison  to  its 
rather  smooth  mode-like  appearance  for  e = 0.4.  The  associated  values 
of  RMS(t')  and  RMS(C')  are  comparable  to  those  for  e=0.4. 

With  the  addition  of  a mean  flow  (y  / 0),  the  FD  model  actually 
becomes  more  accurate  perhaps  reflecting  the  increased  smoothness  of 
the  | field  (Figure  25b).  With  y = 0.  5,  the  FD  model  delivers  a stable 
solution  with  an  accuracy  of  0(20$)  after  5 periods  (Table  I,  case  26). 

For  a mean  flow  of  the  opposite  sense  (y  = -0.  5),  the  errors  are 
comparable  or  slightly  larger.  As  with  y = 0,  the  integrated  errors  grow 
linearly  in  time,  and  the  fields  of  t(r ' and  £'  are  dominated  by  large-scale 
box  mode-like  features.  Similar  remarks  hold  for  Rossby  waves  of 
different  orientation. 
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7.  3 Finite-element  model  results 
(Table  I,  cases  22-37;  Figures  21-23,  25-27) 

The  results  of  the  FE  linear  Rossby  wave  calculations  are  as 
expected  from  Sections  4.  3 and  6.  3.  With  (v,T),  k,  f,  y)  = (32/3.  5,  128, 

3/fn,  2/713,0),  RMS(|')<  RMS(C')  *0(1$)  at  the  end  of  five  periods  of 

integration  (Table  I,  case  22).  In  fact,  RMS(£')  has  nearly  stopped 
increasing  altogether,  though  RMS(if ')  continues  to  rise  quasilinearly. 
Associated  with  these  error  levels  are  fields  of  i|r 1 and  £'  composed 
of  basin- scale  and  grid- scale  features,  respectively. 

With  e >0,  the  FE  model  accumulates  box  mode-like  features  in 
the  streamfunction  field.  Errors  in  the  vorticity  are  resident  at  some- 
what smaller  scales.  Consider,  for  instance,  Figures  21c,  d and  22c,  d 
which  show  the  total  and  perturbation  fields  respectively  at  the  end  of  a 
five -period  integration  with  (v,T) , e,  k,  I,  y)  = (32/3. 5, 64, 0. 4, 3/VT3, 

2/VT3,  0).  For  this  case,  the  FE  model  has  errors  of  0(9$),  a significant 
improvement  over  the  second-order  FD  results  which,  as  remarked, 
have  a large  component  of  spatial  truncation  error.  On  the  contrary,  the 
FE  errors  are  most  sensitive  to  changes  in  r),  at  least  in  the  parametric 
neighborhood  of  our  pivotal  calculation  (Table  I,  cases  23-25). 

Quantitatively  similar  statements  can  be  made  for  simulations  at 
higher  Rossby  number--case  28,  e = 0.8--and  in  the  presence  of  mean 
advection- -cases  26  and  27,  y = + 0.  5.  (Note  that  the  latter  differ  from 
the  nonadvected  Rossby  waves  in  that  they  have  nontrivial  nonlinearities.  ) 

As  with  the  FD  model,  neither  the  increase  in  C nor  the  inclusion  of 
mean  advection  seriously  increases  the  RMS  errors  of  the  FE  model. 

As  a result,  for  constant  v and  q (spat’al  and  temporal  resolution), 

RMS(f')  = 0(1-4$)  and  RMS(£')  = 0(3-10$)  after  5 periods  (Table  I,  cases 
25-28).  The  error  growth  is  quite  consistently  nearly  linear  (Figures  27c,  d) 
with  the  perturbation  streamfunction  appearing  at  the  largest  (basin) 
scales  although  in  a somewhat  less  organized  pattern  than  the  box  mode- 
like features  noted  with  (e,  y)  = (0. 4,  0)--Figures  22d  and  26d. 
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7.4  Pseudospectral  model  results 
(Table  I,  cases  22-37;  Figures  21-27) 

For  0.4s  es  0.  8,  the  spectral  model  suffers  eventual  numerical 
instability  at  some  ts  5 periods.  Figure  24a  shows  a typical  example 
where  (v,T),  e,*Y)  = (16/3.  5,  128,  0.  4,  0).  By  t=1.5  periods,  the  total 
vorticity  is  dominated  by  small-scale  noise;  catastrophic  failure  of  the 
numerical  experiment  occurs  shortly  thereafter.  The  most  intense  grid- 
scale  vorticity  features  occur  at  one  or  more  points  on  the  boundary,  but 
the  noise  is  also  substantial  in  the  interior  along  a line  normal  to  that 
point.  This  is  undoubtedly  due  to  the  nature  of  the  spectral  expansion 
which  ties  points  together  in  just  such  a manner.  The  ultimate  origin  of 
the  PS  instability  is  not  known.  The  site  of  the  instability,  for  instance, 
is  random  and  not  simply  related  to  the  imposed  patterns  of  inflow/outflow 
along  the  domain  margins. 

It  has  been  discovered  empirically,  however,  that  periodic  spectral 
filtering  effectively  controls  the  generation  and  accumulation  of  grid- 
scale  vorticity,  and  prevents  numerical  instability  in  these  nonlinear 
Rossby  wave  experiments.  Figure  24  shows  the  effect  on  one  such  pivotal 
calculation.  By  t = 1.  5 periods,  Q is  entirely  dominated  by  two  regions 
of  high  wavenumber  noise  in  the  unfiltered  calculation.  When  the  simulation 
is  redone,  however,  with  filtering,  the  Rossby  wave  is  easily  advanced  in 
time  to  t = 5 periods.  The  final  field  is  quite  free  of  grid-scale  noise. 

Filtering  of  this  kind  stabilizes  a wide  range  of  nonlinear  Rossby 
wave  calculations  (Table  I,  cases  23-37).  The  resulting  RMS  errors  are 
also  notably  small,  being  no  more  than  a few  percent  for  the  experiments 
recorded  in  Table  I.  The  errors  associated  with  the  filtered  PS  model 
are  typically  many  times  smaller  than  those  of  the  comparable  FD  test 
and  somewhat  smaller  than  those  given  by  the  FE  model.  RMS(^)  and 
RMS(C')  grow  linearly  in  time  (perhaps  with  some  initially  large  value 
of  the  errors  due  to  the  filtering  (Figures  23  and  27)  with  very  little 
accumulation  of  unresolvable  features  in  the  vorticity  field  (Figures  21e,f 
and  25e,f).  The  removal  of  these  small-scale  features  by  filtering  does 
not,  however,  seem  to  have  a strong  effect  on  the  energy  of  the  system. 
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7.  5 Intercomparison 

The  parametric  results  of  the  linear  Rossby  wave  calculations 
confirm  the  conclusions  of  Sections  4.  5 and  6.  5 in  which  it  was  noted 
that  the  FE  and  PS  models  were  in  general  more  accurate  than  the  FD 
code,  except  for  certain  optimal  choices  of  the  computational  parameters 
v and  r).  In  addition,  the  orientation  of  the  Rossby  wave  has  a strong 
influence  on  the  RMS  errors  of  the  FD  model.  All  three  models  are 
stable  for  e = 0. 

When  nonlinearity  is  admitted,  however,  the  finite- difference  and 
finite -element  models  are  alone  capable  of  delivering  stable  and  accurate 
calculations  of  moderate  duration  (t  s 5 periods)  over  a broad  range  of 
parameters.  The  spectral  model  is  typically  unstable  in  these  instances 
unless  it  is  supplemented  by  periodic  spectral  filtering  of  the  vorticity  field. 

As  a preliminary  test  of  the  response  of  the  finite -difference  and 
finite -element  models  to  the  addition  of  a scale  selective  vorticity  filtering 
mechanism,  we  have  redone  the  FD  and  FE  experiments  25,  26  and  28  and 
the  FE  experiments  7 and  21  with  the  application  at  each  time  step  of  a 
16th-order  Shapiro  filter  [33].  The  results  of  these  comparisons  indicate 
that  the  RMS  errors  of  the  FE  calculations  are  generally  lowered  some- 
what by  the  addition  of  filtering  (particularly  RMS(£ ')  whose  smaller 
scale  components  are  being  eliminated  by  the  filtering)  and  its  instabilities 
delayed  (but  not  removed).  The  opposite,  namely  an  increase  of  error 
with  the  application  of  filtering,  is  often  true  of  the  FD  simulations.  It  is 
not  obvious  why  this  should  be  the  case  unless  the  computational  boundary 
condition  used  in  the  FD  formulation  interacts  in  some  systematic  way 
with  the  applied  filtering.  This  possibility  will  be  explored  in  our  next 
series  of  pseudoforecasting  tests  (see  Section  8). 

It  is  of  interest  to  note,  however,  that  all  three  models  (perhaps 
with  some  distribution  of  wavenumber  selective  filtering)  can  be  made  to 
yield  accurate  solutions  to  these  open  domain  problems.  In  fact,  the 
models  have  error  accumulation  characteristics  not  greatly  different  than 
those  noted  in  closed-basin  problems.  Specifically,  the  FE  and  PS  models 
are  many  times  more  accurate  for  given  v and  r)  than  the  FD  model, 
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with  the  PS  being  overall  the  best.  Even  taking  into  account  the 
increased  efficiency  of  the  finite -difference  scheme  (Table  II),  the 
difference  between  the  second-order  and  higher-order  methods  is 
significant.  It  is  estimated  that  the  FE  (N  = 33)  and  PS  (N  = 17)  models 
are,  on  the  average,  15  and  3.  5 times  more  accurate  respectively  than 
a FD  model  with  N ^ 43  for  which  the  running  times  of  all  three  models 
would  be  approximately  equal. 
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8.  CONCLUSION 

We  have  integrated  the  inviscid  barotropic  vorticity  equation  under 
a variety  of  assumed  initial  and  boundary  conditions  corresponding  to 
linear  and  nonlinear  box  modes,  forced  nonlinear  box  modes,  north  wall 
forced  modes  (meander  induced  forcing),  and  linear  and  nonlinear  Rossby 
waves.  The  former  two  classes  of  problems  are  defined  with  closed 
domains;  the  latter  two  are  partially  or  totally  open  with  respect  to  a 
presupposed  external  environment  and  therefore  represent  prototype 
limited-area  calculations  for  the  ocean.  Each  problem  has  been  solved 
using  second-order  finite -difference,  fourth-order  finite -element  and 
infinite-order  spectral  approximation  techniques.  For  each  of  the  three 
models  a series  of  calculations  was  performed  to  determine  its  accuracy, 
stability  and  efficiency  as  a function  of  problem  type  and  the  associated 
physical  and  computational  nondimensional  parameters.  The  most 
important  of  these  parameters  are  e,  the  Rossby  number,  and  v and  r), 
nondimensional  measures  of  the  spatial  and  temporal  resolution  of  the 
numerical  approximation.  The  accuracy  of  model  results  was  determined, 
wherever  possible  by  comparing  to  known  analytic  or  reference  solutions. 
RMS  measures  of  the  errors  in  the  computed  values  of  vorticity,  RMS(C'), 
and  streamfunction,  RMS(f '),  and  a measure  of  the  gain  or  loss  of  globally 
integrated  kinetic  energy,  NDIF  (NRG),  were  tabulated.  Integrations  of 
moderate  length  (5-10  periods  of  the  reference  solution)  were  performed  as 
an  empirical  measure  of  the  functional  dependence  of  model  stability  on 
the  parameters.  As  a result  of  these  calculations,  we  are  able  to  make 
model-model  intercomparative  statements  for  a sequence  of  linear  and 
nonlinear  problems  in  open,  as  well  as  closed,  domains.  To  our 
knowledge,  such  intercomparisons  have  not  previously  been  made.  A more 
lengthy  summary  of  the  parameters  and  error  norms  can  be  found  in 
Section  2.  A complete  discussion  of  results  has  been  given  in  Sections  4-7 
and  Table  I. 

These  tests  have  shown  that  all  three  models  are  capable  of  delivering 
efficient  long-term  solutions  of  acceptable  accuracy  to  linear  and  weakly 
nonlinear  problems  in  both  closed  and  open  domains.  The  results  also 
suggest  that  given  a judicious  selection  of  frictional  (filtering)  mechanism 
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and/or  computational  boundary  condition,  each  of  the  models  can  be  made 
comparably  accurate  for  highly  nonlinear  calculations.  (This  hypothesis 
is  being  tested  in  a related  series  of  experiments).  We  conclude,  there- 
fore, that  any  of  the  physical/numerical  models  investigated  here-- 
modified  perhaps  by  additional  dissipative  or  boundary  condition 
assumptions- -could  be  used  for  the  intended  scientific  applications 
mentioned  in  Section  1. 

Under  the  assumption  of  inviscid  dynamics,  the  operational  per- 
formance of  the  three  models  is  most  sensitively  related  to  the  Rossby 
number,  e.  For  Os  e ^ 0.  2,  all  the  models  are  stable  in  the  long- 
term. Furthermore,  unless  an  optional  choice  of  v and  r|,  the  non- 
dimensional  space  and  time  steps,  is  made,  the  spectral  and  finite- 
element  models  are  the  most  accurate,  and  the  finite- difference  the  least. 
That  this  ranking  reflects  the  formal  spatial  accuracies  of  the  models  has 
been  demonstrated  by  a simple  phase  error  analysis  for  the  linear  box 
mode  problems --Section  4.  2.  The  net  result  of  this  increased  accuracy 
is  that,  for  a given  admis sable  error,  both  the  FE  and  PS  models  are 
many  times  more  efficient  than  the  FD  model  (Section  7.  5).  These 
conclusions  are  valid  independent  of  problem  class. 

Although  the  PS  (and  to  a lesser  degree  the  FE)  models  are 
susceptible  to  eventual  numerical  instability  characterized  by  the 
catastrophic  accumulation  of  grid- scale  vorticity  features,  it  has  been 
found  that  stability  can  often  be  maintained,  and  errors  reduced,  by  a 
periodic  filtering  (smoothing)  procedure.  Consider  for  instance  the  non- 
linear Rossby  wave  experiments  in  the  presence  of  mean  advection.  The 
PS  model  develops  numerical  instabilities  which  appear  as  very  high 
wavenumber  noise  in  the  vorticity  field.  (This  generation  of  noise  may  be 
partially  due  to  our  choice  of  boundary  conditions  - -see  the  following 
remarks.  ) But,  by  selectively  filtering  out  this  grid- scale  vorticity  at 
each  time  step,  the  spectral  model  can  be  made  stable  in  the  long-term 
sense  while  maintaining  a very  high  accuracy  (Table  I,  cases  23-37).  The 
RMS  error  norms  of  the  FE  (but  not  the  FD)  model  are  also  reduced  with 
the  application  of  a scale  selective  vorticity  filter. 
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If  we  had  sufficient  spatial  resolution  (that  is,  scale  separation 
between  the  energy- containing  band  and  our  cutoff  wavenumber,  Kq)  and 
could  devise  the  optimal  filter,  removing  vorticity  in  this  manner  would 
correspond  to  destroying  enstrophy  as  it  reached  K^,  but  to  leaving  energy 
unaffected.  Such  intertial  range  filters  have  been  constructed  on  the  basis 
of  various  turbulence  theory  closure  schemes  [34]  and  are  typically  non- 
linear and  highly  dependent  on  the  energy  spectrum  of  the  field  being 
filtered.  Standard  frictional  mechanisms  such  as  linear  and  second-order 
vorticity  dissipation,  and  the  exponentially  tapered  vorticity  filter  adopted 
here  for  the  spectral  calculations  are  almost  certainly  crude  approximations 
to  these  optimal  filters.  As  we  have  seen,  however,  even  relatively  ad  hoc 
removal  of  small-scale  vorticity  prolongs  the  useful  length  of  integration 
and  increases  the  accuracy  of  highly  nonlinear  simulations. 

Were  we  to  incorporate  a higher  order  frictional  mechanism  into 
our  models  in  order  to  prevent  the  contamination  of  our  results  by  small- 
scale  noise,  however,  the  problem  of  specifying  the  outflow  boundary 
condition  would  take  on  added  importance.  We  already  know  on  the  basis 
of  finite -difference  calculations  detailed  in  Section  4.  2 that  alternate 
specifications  of  can  lead  to  much  smoother  vorticity  fields.  Similar 
findings  for  the  nonlinear  Rossby  waves  (not  presented  here)  also  indicate 
that  certain  outflow  boundary  conditions  are  better  able  to  control  the 
accumulation  of  small-scale  vorticity  near  the  boundaries,  presumably  by 
allowing  the  grid-scale  error  field  to  propagate  more  freely  through  the 
boundaries.  Although  the  optimal  form  of  the  outflow  boundary  condition 
is  a matter  of  some  debate,  there  are  several  strong  candidates  which  have 
proven  useful  in  various  applications.  These  include  the  Sundstronr/Davies 
condition  considered  here  [16,  17],  the  Orlanski  radiation  condition  [11], 
and  various  extrapolatory  techniques  and  multi-dimensional  generalizations 
of  the  method  of  characteristics.  Lastly,  it  is  possible  to  locally  modify 
the  physics  of  the  problem  so  that  small-scale  errors  are  trapped  near  the 
boundaries  and  selectively  dissipated  before  they  can  contaminate  the 
interior  solution  [35,36]. 


A second  series  of  tests  is  now  being  readied  to  investigate  further 
the  effects  of  dissipation  (including  filtering)  and  alternate  specifications 
of  the  outflow  boundary  condition  on  the  accuracy  and  stability  of  related 


limited-area  calculations.  Following  the  methodology  of  this  report,  we 
plan  a systematic  parametric  examination  of  these  effects.  Since  it  is 
expected  that  this  testing  phase  serve  as  a lead-in  to  the  intended 
scientific  applications  outlined  in  Section  1,  the  reference  solution 
selected  for  these  next  tests  will  have  a multiplicity  of  space  and  time 
scales.  Two  approaches  are  possible.  First,  a desired  analytic  solution 
can  be  constructed  (perhaps  assuming  some  known  distribution  of  body 
forces).  Second,  the  open  ocean  problem  can  be  embedded  into  a 
previously  existing  closed-basin  numerical  calculation  of  known  accuracy 
(which  serves  as  the  source  of  boundary  information).  Of  these  two 
approaches,  the  former  is  the  easier  to  implement  but,  in  general,  would 
require  the  assumption  of  some  complicated  time -dependent  forcing 
function.  As  has  been  speculated  in  Section  5,  this  may  significantly 
reduce  the  computational  error  because  the  computed  solution  tends  to  be 
"locked  in"  to  the  forcing.  Using  a forced  analytic  reference  solution  as 
a test  problem  is,  therefore,  neither  physically  nor  numerically  analogous 
to  the  unforced  forecasting/hindcasting  studies  proposed  in  Section  1.  An 
embedding  experiment,  on  the  other  hand,  can  be  made  physically  as  well 
as  numerically  identical  to  a forecasting  study  although  it  requires  that 
boundary  and  verification  data  be  generated  numerically  and  stored  for 
later  access.  For  this  reason,  we  will  pursue  the  embedding  strategy  as 
a basis  for  our  next  series  of  tests. 

The  large-scale  unforced  numerical  simulation  into  which  the  open 
ocean  calculations  will  be  embedded  can  be  conducted  in  either  a closed 
or  periodic  domain.  Solutions  will  be  sought  whose  local  character 
resembles  that  of  the  proposed  scientific  applications;  that  is,  either  the 
evolution  of  a mesoscale  eddy  field  such  as  might  be  observed  in  the  mid- 
ocean or  the  local  dynamics  of  intense  current  regions  might  be  examined 
as  natural  precursors  to  the  MODE/POLYMODE  forecasting  and  EGCM 
jet  instability  studies,  respectively.  Since  it  is  necessary  to  consider 
various  frictional  and  filtering  mechanisms  (including  the  inviscid  limit), 
reference  numerical  solutions  are  needed  for  each  of  these  mechanisms 
although  the  environmental  parameters  could  otherwise  be  held  fixed. 
(Besides  providing  initial  and  boundary  condition  data  for  the  open  ocean 
models,  this  could  provide  an  interesting  study  of  the  effects  of 
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subgrid-scale  parameterization  on  geophysical  turbulence. ) In  addition, 
in  order  to  avoid  biasing  the  results  towards  a particular  numerical 
technique,  the  large-scale  reference  solution  must  be  generated  either 

1)  with  a numerical  model  different  than  those  being  compared,  or 

2)  independently  for  each  numerical  model  being  tested. 

The  pivotal  calculation  for  each  of  these  tests  will  be  inviscid  and 
use  the  Charney- Fjortoft- von  Neumann  closure.  Subsequent  calculations 
will  involve  the  assumption  of  linear  bottom  drag,  high-order  lateral 
friction  and  spectral  filters  (or  their  equivalent  in  physical  space)  in 
various  combinations.  Simultaneously,  we  will  seek  to  minimize  wave 
reflection  from,  and  boundary  layering  near,  open  boundaries  by  invoking 
the  Sundstrom/Davies,  Kreiss,  Orlanski  (or  other)  condition  on  outflow. 

It  will  also  be  of  interest  in  some  cases  to  examine  the  effects  on  the 
accuracy  and  stability  of  the  models  of  adding  errors  of  certain  amplitude 
and  scale  to  the  imposed  boundary  data.  This  will  be  important  in  judging 
the  suitability  of  the  models  for  forecasting  studies  driven  by  observational 
data  in  which  measurement  and  objective  analysis  errors  may  be  large. 

On  this  basis,  it  will  be  possible  to  confidently  select  that  physical  problem/ 
numerical  technique/boundary  condition  set  which  is  most  reliable  for 
limited-area  calculations  with  e»0. 

At  some  point  during  this  sequence  of  tests,  it  will  become 
necessary  to  choose  one  of  the  open  ocean  models  as  that  best  suited  for 
the  intended  scientific  problems.  This  selection  will  be  done  expeditiously 
and  will  be  based  on  such  factors  as  the  efficiency  of  the  models  for  a 
given  accuracy  and  the  generalizability  of  the  models  (to  account  for 
baroclinic  processes  or  alternate  boundary  conditions,  for  instance). 

Once  the  choice  of  numerical  technique  has  been  made,  a baroclinic 
(two- level)  version  of  barotropic  model  will  be  constructed  and  readied  for 
use.  Some  limited  testing,  along  the  lines  of  this  report,  will  be 
necessary  to  demonstrate  that  our  conclusions  based  on  these  barotropic 
calculations  do  indeed  carry  over  to  the  baroclinic  problem. 
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APPENDIX  I 

Implementation  of  the  Sundstrom/Davies  Boundary  Condition 

Consider  a region  near  the  boundary  with  the  following  local 
ordering. 
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The  vorticity  equation  applied  at  the  point  0 will  involve  vorticity 
values  on  the  other  eight  numbered  points. 
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Where  Q = £At/h  and  Jr  is  the  Arakawa  Jacobian  term  for  the  r-th 
point  evaluated  at  point  0. 

The  Sundstrom/Davies  closure  for  the  boundary  point  2 is 


rk+l  k- 1 _ k k 

^>0  ~ *2  ^6 


(A2) 


k+1 

can  be  eliminated  from  (Al)  and  (A2)  and  the  result  can  be  written: 

(A3) 


QJ1£$:+  (1  + QJ2)Ck  + QJ3C3  = R2 


where 
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In  this  numerical  procedure,  the  interior  vorticity  and  streamf unction 
are  computed  before  the  boundary  vorticity.  Equation  (A3)  is  arranged  so 


that  contains  only  known  quantities.  For  each  point  along  a boundary 

(but  not  at  corners),  (A3)  applies. 


Were  it  not  for  corner  points,  the  system  of  equation  (A3)  would  be 
of  tridiagonal  form  and  easily  invertable.  Consider  a region  near  the 
(southeast)  corner. 
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For  the  boundary  condition  at  the  corner,  we  take  the  spatial 
average  along  the  diagonal. 


.k+1  . rk-l* 
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It  would  seem  that  near  corners  a pentadiagonal  system  is  required  since 
the  Jacobian  evaluated  at  0 involves  the  five  unknown  boundary  values  at 
points  one  through  five.  However,  the  Sundstrom/Davies  conditions  at 
points  two  and  four  are 

rk+1  + rk_1  - rk  + rk  ,An 

*»0  + ^0  “ ^2  + ^6  (A6) 

and 
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Thus, 
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and 
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This  means  that  boundary  points  which  are  corner  neighbors  can  be 
expressed  in  terms  of  the  corner  vorticity  and  known  values.  (A8)  and 
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(A9)  are  applied  in  conjunction  with  the  equations  requiring  vorticity 
values  of  the  corner  neighbor,  i.  e.  , the  corner  and  the  two  boundary 
points  2h  from  the  corner. 

The  unknowns  in  the  tridiagonal  system  are  given  a cyclic  ordering 
excluding  the  corner  neighbor  points. 


In  the  diagram  the  ordering  is  started  in  the  southwest  corner,  and 
(A3)  is  a tridiagonal  system  with  cyclic  boundary  conditions  since  point  1 
and  4K  are  connected.  By  a reordering  of  the  points,  the  tridiagonal 
cyclic  system  can  be  transformed  into  a pentadiagonal  system  which 
requires  twice  the  computation  time  of  a tridiagonal  system. 

At  inflow  points,  the  vorticity  is  known  and  (A3)  is  replaced  by 

= Cg  (A10) 

where  Cg  is  the  specified  value.  If  the  origin  of  the  ordering  (the  south- 
west corner  in  the  diagram)  were  inflow,  the  off- tridiagonal  terms 
expressing  the  cyclic  nature  vanish  and  we  are  left  with  a simple 
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tridiagonal  system.  If  the  southwest  corner  is  not  inflow,  renumber  the 
boundary  point 

i'  = (i  + L - 1)  mod(4K)  + 1 

where  L is  an  inflow  point.  Finally,  the  values  at  the  eight  corner 
neighbor  points  are  obtained  from  (A8)  and  (A9),  or  specified  if  inflow. 
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TABLE  CAPTIONS 

Table  I:  Maximum  RMS  error  measures --RMS(^i '),  RMS(Q'),  and 

NDIF  (NRG)- -for  the  finite -difference  (FD),  finite-element 
(FE)  and  pseudospectral  (PS)  models  as  a function  of  problem 
class,  duration  of  experiment,  and  the  associated  nondimen- 
sional  parameters.  The  first  seven  columns  refer  to  the 
experiment  number  and  the  quantities  e,  A,  x^,  N,  v and  r| 
which  are,  respectively,  the  Rossby  number,  the  number  of 
half  wavelengths  of  the  reference  solution  within  the  domain, 
the  nondimensional  basin  size,  the  number  of  spatial  degrees 
of  freedom  in  each  direction  and  nondimensional  measures  of 
the  spatial  and  temporal  resolution.  Intermediate  columns 
are  reserved  for  problem- dependent  parameters  which  have 
been  introduced  in  Sections  4-7. 

Table  II:  Approximate  model  running  times  as  a function  of  N (CPU 

time  in  seconds  on  the  NCAR  CDC  7600  per  100  time  steps). 
The  ratios  of  the  running  times  (also  listed)  indicate  that  the 
computational  time  increases  as  approximately  N for  the 
FD  and  FE  models  and  as  Nin  N for  the  PS  model. 
Running  times  for  the  other  linear  and  nonlinear  model 
problems  are  comparable  to  those  quoted  here  for  the  Rossby 
wave  calculations. 


Linear  and  Nonlinear  Box  Modes 


Exp't  £ A X N v n m n Duration  RMS(^')  RMS(C')  NDIF  (NRG)  Comments 


0.86  -57.  x 10  filter 


Exp't  £ A X_  N V n case  Duration  RMS  W)  RMS  {£')  NDIF  (NRG)  Comments 


Exp't  e A X N V n k #.  Y Duration  RMS(^')  RMS(C')  NDIF  (NRG)  Comments 


Exp’t  £ A X N V nk  a Y Duration  RMS(\Jj')  RMS ( £ ' ) NDIF  (NRG)  Comments 


TABLE  II 


Linear  Rossby  Waves 


Nonlinear  Rossby  Waves 


N 

FD 

FE 

PS 

FD 

FE 

PS 

17 

1.5 

3.1 

3.7 

2.0 

3.4 

12.1 

33 

6.3 

13.1 

10.6 

8.3 

14.0 

31.1 

33/17 

4.2 

4.2 

2.9 

4.2 

4.1 

2.6 

Note  that  contour  plots  and  RMS  error  curves  are  not  necessarily  scaled 
similarly,  even  within  a single  figure. 

Figure  1:  Linear  box  mode  solution  at  t=10  periods  with  (r|,m,n)  = (128,1,1). 
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Figure  2:  as  Figure  1 

(a,b)  FD  model,  v = 16^?,  Cl  = (0.08,  0.02) 

(c,d)  (^,,^’):  FE  model,  v = 16\/2,  Cl  = (4.0,  1.0)  x 10_3 

(e,f)  (£', !/>'):  PS  model,  v=  8\^,  Cl  = (8.0,  3.0)  x 10*3 


Figure  3:  as  Figure  1 

(a,b)  RMS ( £ ' ,\p' ) : FD  model,  v = 16 /?»  0 < t < 10  periods 

(c,d)  RMS(^‘ ,t|)')  s FE  model,  v = 16^,  0 < t < 5 periods 

(e,f)  RMS ( £ * ,ip')  : PS  model,  v = 8^J,  0 < t < 10  periods 


Figure  4:  Nonlinear  box  mode  solution  at  t=5  periods  with  (v»lVm'n'£)  = 

(16\/7,  64,  1,  1,  0.2). 

(a,b)  (£,^):  FD  model.  Cl  = (0.3,  0.1) 

(c,d)  (^,rp)  : FE  model.  Cl  = (0.3,  0.1) 

(e , f ) (£,\p)  : PS  model,  Cl  = (0.4,  0.1) 


Figure  5: 


Figure  6: 


as 

Figure 

4 

(a. 

b) 

<5* 

: FD 

(c. 

d) 

<C’ 

rip') 

: FE 

(e, 

f) 

<c* 

rlp’) 

: PS 

as 

Figure 

4, 

Cl  = 

(a) 

FD 

model 

(b) 

c'* 

FD 

model 

(c) 

FD 

model 

Cl  = (0.1,  0.03) 

8.0  x 10-3) 


Cl  = (0.1, 

Cl  = (0.4,  0.03) 


Se  = 5a 


specified  by  Kreiss  condition 


£ specified  by  Sundstrom/Davies  condition 

u 


J 
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Figure  7:  as  Figure  4 


(a,b) 

(C*4»)  * 

PS 

model, 

t=2 

periods. 

Cl  = (0. 

(c,d) 

(C^)s 

PS 

model. 

t=2 

periods. 

filtered 

(e,f ) 

PS 

model. 

t=5 

periods. 

filtered 

0.1) 


0.1) 

0.1) 


Figure  8:  Forced  nonlinear  box  inodes  with  (v»n»a»b»e)  = (32,128,l/v/7,l/\/7,0.2) 


(a,b) 

(£,!/>)  : 

FD  model. 

t=2 

periods. 

Cl  = (0.4, 

0.1) 

(c,d) 

<?»*> : 

FE  model. 

t=5 

periods. 

Cl  = (0.7, 

0.1) 

(e,  f ) 

<?.<!»>  » 

PS  model. 

t=5 

periods. 

Cl  = (0.4, 

0.1) 

Figure  9:  as  Figure  8 


(a,b) 

FD 

model , 

t=2 

periods, 

Cl  = (0.04,  1.0 

x 10 

(c,d) 

<C\<K>  : 

FE 

model , 

t=5 

periods. 

Cl  = (0.1,  4.0 

x 10 

(e,f) 

(?’, (I/'): 

PS 

model. 

t=5 

periods. 

Cl  = (0.06,  2.0 

x 10 

3 


Figure  10:  as  Figure  8 

(a,b)  RMS(c',^'):  FD  model,  0 < t < 2 

(c,d)  RMS(c',ip'):  FE  model,  0 < t < 5 

(e,  f ) RMS ( £ ' ,\j> ' ) : pS  model,  0 < t < 5 


Figure  11:  Case 

1 linear 

north  wall 

forced  modes  at 

(a,b) 

FD  model. 

Cl  = (0.3,  0.1) 

(c,d) 

FE  model, 

Cl  = (0.3,  0.1) 

(e,f) 

PS  model. 

Cl  = (0.3,  0.1) 

periods 

periods 

periods 


(V.n)  = (32/3,  64) 


Figure  12:  as  Figure  11. 


(a,b) 

(C'.l|»,>: 

FD  model. 

Cl  = (0.1, 

(c,d) 

(C',lK): 

FE  model. 

Cl  = (0.01 

(e,f) 

(c‘,li»')  : 

PS  model, 

Cl  = (0.01 

Figure  13:  as  Figure  11, 

0 < t < 5 

periods 

8. 0x  lO-3) 
8. Ox  IQ-3) 


(a,b)  RMS  ( £ 1 ,ijj' ) 
(c,d)  RMS(^',^') 
(e,f)  RMS(£\4,') 


FD  model 
FE  model 
PS  model 


i 
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Figure  14:  Case  2 linear  north  wall  forced  modes  at  t=5  periods,  with 

(V,n)  = (64/3,  64) 


(a,b) 

FD  model, 

Cl  = (10.0, 

3.0) 

(c,d) 

(Z.ifi): 

FE  model, 

Cl  = (10.0, 

3.0) 

(e,  f) 

(£»♦>  : 

PS  model, 

Cl  = (10.0, 

2.0) 

Figure  15:  as  Figure  14 


(a,b) 

(C,,<f/’>  : 

FD 

model. 

Cl  = 

(2.0, 

0.06) 

(c,d) 

(C’^’)  : 

FE 

model , 

Cl  = 

(0.5, 

0.01) 

(e,f ) 

(SW): 

PS 

model , 

Cl  = 

(0.1, 

0.01) 

Figure  16:  as  Figure  14,  0 < t < 5 periods 

(a,b)  RMS(£',^’):  FD  model 
(c,d)  RMS  ( £ 1 , 1 ) : FE  model 
(e , f ) RMS  (£ ' ,ijj ' ) : PS  model 


Figure  17:  Case  1 nonlinear  north  wall  forced  modes  at  t=5  periods  with 

(V,rj,e)  = (32/3,64,0.2) 


(a,b) 

FD  model. 

Cl  = (2.0, 

0.1) 

(c,d) 

(?,^): 

FE  model , 

Cl  = (2.0, 

0.1) 

(e,f) 

(t;,i|>) : 

PS  model. 

Cl  = (2.0, 

0.1) 

Figure  18:  as  Figure  17,  t=0.5  periods 


(a,b)  u\r) 
(c , d)  (£'  ,\J/') 
(e,f ) (S',r) 

FD  model.  Cl  = 
FE  model,  Cl  = 
PS  model.  Cl  = 

(0.06,  0.01) 

(0.04,  0.01) 

(0.04,  0.01) 

Figure  19: 

as  Figure  17, 

t=5.0  periods 

(a,b)  (C',**) 
(c,d)  (C',^') 
(e,f)  (C.ifi') 

FD  model,  Cl  = 
FE  model.  Cl  = 
PS  model.  Cl  = 

(2.0,  0.1) 

(2.0,  0.1) 

(2.0,  0.07) 

Figure  20: 

as  Figure  17 

(a,b)  (?,!/;): 
(c,d)  (^,lf»): 
(e,f ) (C,^): 

PS  model,  t=l . 5 
PS  model,  t=l . 5 
PS  model,  t=5.0 

periods,  Cl  = (3.0,  0.1) 
periods,  filtered,  Cl  = 
periods,  filtered.  Cl  = 

(0.9,  0.1) 

(2.0,  0.1) 
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Figure  21:  Nonlinear  Rossby  waves  at  t=5  periods  with  (v,Ti*k»^/e»y)  = 

(32/3.5,  64,  3/ yd 3,  2/VTT,  0.4,  0.0) 

(a,b)  (C,l^):  FD  model,  Cl  = (0.2,  0.2) 

(c,d)  (£,i jj)  : FE  model.  Cl  = (0.2,  0.2) 

(e,f)  ( ^ ,\J>)  s PS  model,  CX  = (0.2,  0.2),  filtered 


Figure  22:  as  Figure  21 


(a,b) 

(C,  «'>* 

FD  model, 

Cl  = (0.04,  0.04) 

(c,d) 

(C'.^') * 

FE  model. 

Cl  = (0.05,  0.02) 

(e,f) 

PS  model, 

Cl  = (8.0x10-3,  6.0  xlO-3),  filtered. 

Figure  23:  as  Figure  21,  0 < t < 5 periods 

(a,b)  RMS(C',ifj')  : FD  model 
(c,d)  RMS  (£  1 ,lj/ ' ) : FE  model 
(e,f)  RMS (£ ' ,iji  * ) : PS  model,  filtered 


Figure  24:  as  Figure  21 

(a,b)  (£,ij;):  PS  model,  t=2.5  periods.  Cl  = (0.6,  0.2) 

(c,d)  (?,li»)  : PS  model,  t=2.5  periods,  filtered.  Cl  = (0.2,  0.2) 

(e,f ) (£,\p):  PS  model,  t=5.0  periods,  filtered.  Cl  = (0.2,  0.2) 


Figure  25:  Nonlinear  Rossby  waves  with  (V,rukf Jl^y)  = (32/3.5,64,3/\/T3,2/\/T3, 

0.4,0. 5) 


(a,b) 

U,<iO  : 

FD 

model , 

t=5 

periods. 

Cl  = (0.2, 

0.7) 

(c,d) 

(5rW  * 

FE 

model , 

t=5 

periods, 

Cl  = (0.2, 

0.7) 

(e,f ) 

PS 

model , 

t=5 

periods. 

Cl  = (0.2, 

0.7) , filtered 

Figure  26:  as  Figure  25 


(a,b) 

<c’,r> : 

FD  model, 

t=5 

periods. 

Cl  = (0.07, 

0.03) 

(c,d) 

FE  model, 

t=5 

periods , 

Cl  = (0.04, 

0.01) 

(e,f) 

PS  model. 

t=5 

periods. 

Cl  = (0.01, 

5. Ox  10 

filtered 


Figure  27:  as  Figure  25,  0 < t < 5 periods 

(a,b)  RMS(C  ,</>')  : FD  model 
(c,d)  RMS ( £ ' ,\j)' ) : FE  model 
(e,f)  RMS(£',iK):  PS  model,  filtered 
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